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1.  Executive  Summary 


a.  Phase  I  Program  Objectives1 


1 .  Demonstrate  the  feasibility  of  formulating  a  nontraditional  model  of  LRE  stability 
in  terms  of  well-known  thermomechanical  system  principles,  with  the  potential  to 
provide  quantitative  predictions  of  combustion  stability  in  an  idealized  engine 
configuration. 

2.  Demonstrate  that  the  thermomechanical  approach  can  be  used  to  establish  a 
direct  correspondence  between  distributed,  transient  heat  addition  and 
mechanical  (wave)  disturbance  growth  in  a  semi-confined  geometry. 

3.  Formulate  stability  concepts  attributable  to  thermomechanics  in  a  semi-confined 
system. 


During  the  nine  month  Phase  I  contract  period  (15June2010-14March201 1)  the  research 
program  has  focused  on: 

•  an  evaluation  of  the  flamelet  concept  application  as  currently  implemented  to 

supercritical  turbulent  flows  [1],  and 

•  a  quantitative  thermomechanical  formulation  for  gas  motion  and  acoustic 

disturbances  induced  by  localized,  spatially  distributed,  transient  heat 
addition,  including  that  from  exothermic  chemical  reactions,  as  a  metaphor 
for  physical  and  chemical  phenomena  occurring  in  a  liquid  rocket  engine 
(LRE)  combustion  chamber, 

to  demonstrate  the  feasibility  of  the  objectives  cited  above.  Summaries  of  the  results 
obtained  and  their  technical  value  to  the  AFRL  are  given  in  b.  below,  based  on  the  detailed 
technical  papers  in  the  Appendices. 


1  From  the  Phase  I  proposal  to  AFOSR  (Sept.  2009) 
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b.  Phase  I  Research  Accomplishments 
b.l-Flamelet  concept  application  to  supercritical  flows 

The  focus  of  this  study  was  on  the  examination  of  the  subgrid-scale  (SGS)  scalar 
variance,  as  related  to  the  potential  utilization  of  fiamelet  models  under  supercritical 
conditions. 

We  first  addressed  the  issue  of  the  SGS  scalar  variance,  oz ,  as  related  to  the  LES 
modeling  of  the  turbulent  reaction  through  the  fiamelet  model.  Although  we  do  not  have 
yet  a  database  of  reactive  turbulent  supercritical  flows,  since  the  fiamelet  model  relies 
on  the  definition  of  a  conserved  scalar  Z,  we  have  the  necessary  information  from  our 
previous  species-mixing  study  to  perform  the  investigation;  if  a  conserved  scalar  can  be 
defined  for  the  reactive  flow  situation,  because  of  the  additional  coupling  due  to 
reaction,  this  scalar  is  expected  to  have  a  more  complex  behavior  than  the  present  one 
resulting  only  from  mixing,  so  the  results  of  this  study  are  conservative  in  terms  of  the 
level  of  difficulty  to  model  the  scalar  variance. 

In  a  first  step  of  the  study,  we  derived  the  oz  equation  [1,  2]  in  two  forms.  The  first 
form  highlighted  the  new  terms  in  the  equation  with  respect  to  the  atmospheric-pressure 
conditions  case,  and  the  second  form  emphasized  the  new  SGS  terms  in  the  equation 
that  would  require  modeling,  as  well  as  identified  the  nature  of  each  term  (i.e. 
production  or  dissipation).  The  result  of  assessing  the  activity  of  all  terms  in  the  two 
equations  was  to  identify  two  new  SGS  terms  that  are  have  comparable  magnitude  to 
the  classical  terms  and  thus  cannot  be  neglected:  (i)  a  SGS  diffusion  term,  and  (ii)  a 
SGS  Soret  term  that  has  a  dissipative  effect.  Given  the  lack  of  models  for  these  terms,  it 
was  decided  that  trying  to  model  them  in  the  oz  equation  for  the  purpose  of  solving  the 
equation  was  bound  to  be  an  unrewarding  approach,  and  instead,  it  was  decided  to 
propose  direct  models  for  ctz.  Examination  of  the  database  over  several  binary-species 
systems  revealed  that  indeed  the  PDF  of  Z  is  best  approximated  by  a  /3  PDF,  which 
justified  investing  time  in  modeling  oz  since  the  PDF  can  be  constructed  from  its  mean 
provided  by  the  LES  solution  and  a  model  of  oz. 

Thus,  in  the  second  step  of  the  study  we  explored  two  models  for  oz\  (1)  the 
approximate  deconvolution  model  (ADM),  and  (2)  a  dynamic  gradient-based  (GRD) 
model.  These  models  are  conceptually  different  as  ADM  relies  on  a  mathematical 
derivation  with  no  assumptions  required  regarding  the  nature  of  the  scale  interactions, 
while  the  GRD  relies  on  the  mixing-length  hypothesis  in  conjunction  with  an  equilibrium 
assumption.  The  ADM  is  generally  very  good  to  excellent  at  relatively  small  filter  widths, 
but  as  the  filter  width  becomes  large,  there  is  evidence  that  the  series  expansion  on 
which  ADM  relies  does  not  converge  with  increasing  order  of  approximation,  as  an 
overshoot  becomes  increasingly  visible  with  increasing  filter  width.  Thus,  ADM  should 
be  considered  here  as  an  approximation  rather  than  a  converging  procedure  to  obtain 
oz,  and  in  this  particular  case  an  approximation  order  of  3  is  recommended.  For  the 
GRD-based  model,  we  formulated  a  new  model  which  we  exercised  and  compared  to 
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the  results  obtained  with  a  classical  model  in  conjunction  with  two  approximations  for 
the  filter  width.  Compared  to  the  classical  model  which  increasingly  departs  from  the 
filtered  DNS  template  with  either  type  of  filter  approximation,  our  model  maintains  high 
fidelity  even  at  very  large  filter  widths  and  in  a  range  where  test  filtering  is  performed 
close  to  the  production  range  of  wavenumbers. 

A  manuscript  has  been  submitted  for  publication  [1]  and  a  presentation  [2]  has  been 
made  at  a  national  conference  of  this  subject. 

b.2-Thermomechanics  of  reactive  gases 

Transient,  spatially  distributed  combustion  in  a  turbulent  flow  is  thought  to  be  the 
source  of  acoustic  and  stronger  mechanical  disturbances  in  a  LRE  chamber. 
Surprisingly,  a  quantitative  model  for  mechanical  wave  generation  in  a  transient, 
spatially  distributed  reacting  flow  does  not  appear  to  be  available  in  the  technical 
literature.  As  a  first  step  in  the  development  of  such  a  theory,  Kassoy  [3]  describes  a 
thermomechanics-based  study  to  identify  a  quantitative  relationship  between  time 
resolved,  localized  power  deposition  into  a  finite  volume  of  inert  gas  and  the  magnitude 
of  acoustic  or  stronger  disturbances  generated  in  the  neighboring  unheated,  unconfined 
gas.  The  results  demonstrate  that  the  thermomechanical  response  of  the  heated  gas 
depends  on  a  complex  relationship  between  the  magnitude  of  energy  deposited  and  the 
time  scale  for  that  deposition,  but  not  simply  on  the  power  deposition  alone.  The  theory 
is  limited  to  heating  time  scales  small  compared  to  the  acoustic  time  of  the  volume,  a 
condition  necessary  for  the  pressure  to  rise  with  temperature  (while  the  density  change 
is  marginal).  When  the  energy  deposition  is  quantitatively  limited,  nearly  constant 
volume  heating  occurs  (near-inertial  confinement),  characterized  by  a  small  internal  gas 
expansion  Mach  number,  defined  with  respect  to  the  high  temperature  hot  spot  speed  of 
sound.  Temperature  and  pressure  rise  together  within  the  hot  spot.  Localized  high 
pressure,  relative  to  that  in  the  cold  neighboring  gas,  is  necessary  to  create  the 
pressure  gradient  source  of  induced  mechanical  motion  (both  fluid  speed  and  acoustic 
disturbances).  Gas  expelled  from  the  boundary  of  the  hot  spot  (the  “piston  effect”)  is 
the  source  of  mechanical  disturbances  in  the  unheated  environmental  gas.  The 
expelled  gas  Mach  number,  measured  with  respect  to  the  cold  gas  speed  of  sound  will 
be  exceptionally  small  when  the  energy  deposition  is  sufficiently  small  relative  to  the 
quantitative  limitation  referred  to  above.  The  resulting  mechanical  disturbances  in  the 
cold  environment  are  linear  acoustic  waves.  Larger  energy  deposition  within  the 
quantitative  limitation  is  associated  with  much  larger  expelled  gas  Mach  numbers  and 
significantly  stronger  shock  wave  propagation  in  the  cold  environmental  gas.  Beyond 
the  quantitative  limitation  referred  to  above,  the  heating  process  is  fully  compressible, 
characterized  by  a  very  large  internal  Mach  number.  The  asymptotic  theory 
demonstrates  that  this  type  of  thermomechanical  response  to  energy  deposition  is  the 
source  of  strong  blast  waves  [4]  in  the  unheated  external  environment.  The  formulation 
in  Ref.  3  has  been  generalized  during  the  Phase  I  program  to  examine 
thermomechanical  effects  in  a  reactive  gas  with  one-step  kinetics  [5].  The  complete 
technical  description  appears  in  Appendix  3.b 
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The  primary  objective  of  the  study  in  Ref.  5  is  to  quantify  the  thermomechanical 
response  of  a  reactive  gas  affected  by  a  localized  exothermic  chemical  reaction. 
Relative  to  the  objectives  of  the  Phase  I  program  cited  above,  this  type  of  theory 
provides  fundamental  insights  into  the  sources  of  mechanical  disturbances  within  a  LRE 
chamber  containing  a  turbulent  reacting  flow.  Crucial  length  and  time  scales,  as  well  as 
dominant  physico-chemical  processes  derived  from  an  asymptotic  formulation  may  be 
employed  by  numerical  modelers  to  fully  resolve  transient  reactive  flow  phenomena  and 
to  produce  reliable  predictions  of  combustion  stability  in  realistic  LRE  configurations. 

The  theory  in  Ref.  5  is  formulated  for  a  subcritical  gas  with  one-step,  high  activation 
energy  exothermic  kinetics.  The  asymptotic  formulation  demonstrates  that  the 
thermomechanical  response  of  a  reactive  gas  to  localized,  rapid2  chemical  heat  addition 
is  similar  in  nature  to  that  for  an  inert  gas  with  external  heating  [3].  Nearly  inertially- 
confined  heating  occurs  only  when  the  chemical  heat  release  and  the  high  activation 
energy  are  quantitatively  restricted.  Within  that  limitation  one  finds  linear  acoustic  wave 
generation  for  the  smallest  range  of  energy  addition  and  stronger  wave  generation  for 
more  significant  energy  deposition. 

Numerical  solutions  to  the  equations  describing  the  nearly-inertially  confined 
reaction-generated  heat  addition  process  show  that  a  spatially  distributed  reaction  wave 
appears  spontaneously  in  the  hottest  portion  of  the  hot  spot,  and  propagates  through 
the  relatively  slowly  moving  fluid  at  a  supersonic  speed,  relative  to  the  hot  gas  speed  of 
sound.  As  the  front  nears  the  much  colder  boundary  of  the  hot  spot  it  decelerates 
significantly  and  steepens  considerably.  The  configuration  evolves  toward  a 
discontinuous  front  separating  hot,  high-pressure,  burned  gas  on  one  side  from  cold, 
low  pressure  reactant  on  the  other  side.  Although  the  pressure,  temperature  and 
reactant  concentration  jumps  across  the  front  are  reminiscent  of  a  reactive  shock  wave, 
the  front  does  not  propagate  at  a  supersonic  speed  relative  to  the  unburned  cold 
reactant  near  the  boundary  of  the  hot  spot.  In  fact,  during  the  relevant  heating  time 
scale  the  entire  process  occurs  in  a  nearly  incompressible  medium.  Fluid  expelled  from 
the  hot  spot  boundary  acts  as  source  of  mechanical  disturbances  propagating  into  the 
neighboring  cold  gas.  The  amplitude  of  those  disturbances  depends  upon  the  energy 
addition  level  during  the  reactive  phase  of  the  hot  spot. 

In  analogy  with  the  inert  gas  results,  a  fully  compressible  heat  addition  process 
occurs  when  the  chemical  energy  addition  exceeds  the  quantitative  restriction  referred 
to  above.  The  internal  expansion  Mach  number,  defined  relative  to  the  hot  gas  speed  of 
sound  is  substantial  and  can  be  supersonic.  The  hot  spot  transients  must  be  determined 
from  a  numerical  solution  to  the  compressible,  reactive,  gasdynamic  equations,  to  be 
completed  in  the  future. 

In  summary,  an  asymptotic  formulation  of  the  reactive  hot  spot  problem  has  been 
developed  to  quantify  the  thermomechanical  response  of  reactive  gas  to  localized 
energy  addition  on  a  time  scale  short  compared  to  the  acoustic  time  of  the  finite  volume. 


2  Heating  time-scale  short  compared  to  the  local  acoustic  time 
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Explicit  results  have  been  obtained  for  relatively  modest  levels  of  energy  addition 
compatible  with  nearly-constant  volume  heating  and  linear  acoustic  wave  propagation  in 
the  gas  beyond  the  hot  spot.  The  formulation  includes  a  rational  derivation  of  the 
compressible,  reactive,  gasdynamic  equations  describing  the  thermomechanical 
response  of  a  hot  spot  to  larger  values  of  energy  addition.  This  comprehensive 
quantitative  analysis  substantiates  the  feasibility  of  using  thermomechanical  concepts  to 
achieve  many  of  the  objectives  cited  in  la  above. 
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Abstract 

To  model  the  Subgrid-Scale  (SGS)  scalar  variance  under  supercritical-pressure  conditions,  an 
equation  is  first  derived  for  it.  This  equation  is  considerably  more  complex  than  its  equivalent  for 
atmospheric-pressure  conditions.  Using  a  previously  created  Direct  Numerical  Simulation  (DNS) 
database  of  transitional  states  obtained  for  binary-species  systems  in  the  context  of  temporal  mix¬ 
ing  layers,  the  activity  of  terms  in  this  equation  is  evaluated  and  it  is  found  that  some  of  these 
new  terms  have  magnitude  comparable  to  that  of  governing  terms  in  the  classical  equation.  Most 
prominent  among  these  new  terms  are  those  expressing  the  variation  of  diffusivity  with  thermo¬ 
dynamic  variables  and  a  Soret  terms  having  dissipative  effects.  Since  models  are  not  available  for 
these  new  terms  that  would  enable  solving  the  SGS  scalar  variance  equation,  the  adopted  strategy 
is  to  directly  model  the  SGS  scalar  variance.  Two  models  are  investigated  for  this  quantity,  both 
developed  in  the  context  of  compressible  flows.  The  first  one  is  based  on  an  approximate  deconvo¬ 
lution  approach  and  the  second  one  is  a  gradient-like  model  which  relies  on  a  dynamic  procedure 
using  the  Leonard  term  expansion.  Both  models  are  successful  in  reproducing  the  SGS  scalar 
variance  extracted  from  the  filtered  DNS  database,  and  moreover,  when  used  in  the  framework  of 
a  Probability  Density  Function  (PDF)  approach  in  conjunction  with  the  /3-PDF,  they  excellently 
reproduce  a  filtered  quantity  which  is  a  function  of  the  scalar.  For  the  dynamic  model,  the  pro¬ 
portionality  coefficient  spans  a  small  range  of  values  through  the  layer  cross-stream  coordinate, 
boding  well  for  the  stability  of  Large  Eddy  Simulations  using  this  model. 
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I.  INTRODUCTION 


One  of  the  most  useful  concepts  in  turbulent  combustion  is  that  of  a  ‘flamelet’  ([1]).  Its 
usefulness  stems  from  the  fact  that  the  flamelet  model  is  applicable  in  conjunction  with  finite 
rate  chemistry  and  it  is  the  only  turbulent  combustion  model  which  contains  the  essential 
coupling  between  classically-computed  molecular  mass- diffusion  and  chemical  reaction  rate. 
This  model  is  applicable  when  the  chemical  reactions  have  a  relatively  short  characteristic 
time  with  respect  to  the  flow,  so  that  the  flame  is  thin  and  turbulent  eddies  move  it  and 
distort  it  as  an  entity  rather  then  penetrating  into  the  flame  and  affecting  there  the  reaction. 
When  the  flame  is  not  penetrated  by  turbulence,  the  interior  of  the  flame  is  akin  to  a  laminar 
flame,  and  turbulence  modeling  is  only  needed  for  the  flow  external  to  the  flame;  this  is  a 
substantial  advantage  for  modeling  of  turbulent  reactive  flows.  Finding  the  evolution  of  the 
flame  is  then  reduced  to  obtaining  its  statistical  position.  For  non-premixed  combustion, 
several  assumptions  and  mathematical  manipulations  of  the  governing  equations  lead  to  an 
equation  devoid  of  chemical  sources  for  a  quantity  which  is  thus  a  conserved  scalar.  This 
quantity  is  called  ‘the  mixture  fraction’,  Z.  and  the  statistics  of  the  stoichiometric  value 
of  T,  Zs.  determine  the  flame  location;  whereas  the  temporal  solution  of  the  species  mass 
fraction,  Y.  and  temperature,  T,  in  the  Z  coordinate  system  determines  the  internal  flame 
structure.  In  the  species  and  energy  equations  written  in  the  Z  coordinate  system,  the 
scalar  dissipation,  x,  acts  as  a  diffusion  coefficient.  This  scalar  dissipation  is  that  within 
the  flamelet,  and  thus  is  at  a  scale  smaller  than  the  solution  of  the  flow  field;  it  is  thus  the 
Subgrid-Scale  (SGS)  scalar  dissipation.  Generally,  it  is  hypothesized  that  Z  assumes  the 
form  of  a  3  Probability  Density  Function  (PDF)  ([2,  3]),  and  thus  finding  this  distribution 
reduces  to  obtaining  its  mean  and  variance.  Summarizing,  in  order  to  utilize  the  flamelet 
model  one  must  have  available  the  mean  and  variance  of  Z.  and  the  SGS  y.  In  the  context 
of  Large  Eddy  Simulation  (LES)  in  which  the  filtered  governing  equations  are  solved  subject 
to  models  included  to  introduce  palliatives  for  the  filtered-out  SGS  effects,  the  mean  Z 
value  is  known  from  the  LES  solution.  However,  neither  the  SGS  Z  variance  nor  the  SGS 
dissipation  of  Z  can  be  found  from  the  LES  solution,  and  thus  they  must  be  modeled.  For 
example,  several  investigators  [4,  5]  modeled  the  SGS  scalar  variance  as  proportional  to 
the  resolved  (i.e.  LES)  scalar  dissipation,  a  model  which  follows  from  the  classical  mixing 
length  assumption  considering  that  the  scalar  variance  is  a  measure  of  the  level  of  scalar 
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fluctuations. 

Given  the  aforementioned  importance  of  the  SGS  scalar  variance  in  modeling  combustion, 
and  the  numerous  studies  devoted  to  solving  the  SGS  scalar  variance  equation  by  modeling 
those  of  its  terms  that  are  not  computable  from  the  LES  solution  (e.g.  [6,  7]),  we  are  here 
investigating  the  modeling  of  the  SGS  scalar  variance  under  supercritical  pressure  conditions 
and  use  a  previously  created  Direct  Numerical  simulations  (DNS)  database  for  examining 
the  capability  of  several  models.  The  difference  between  this  situation  and  the  much  studied 
atmospheric-pressure  conditions  is  that  here  the  species  mass  fluxes  include  Soret  effects, 
that  real  gas  equations  of  state  (EOSs)  instead  of  the  perfect  gas  EOS  are  used  and  that  the 
transport  properties  are  all  complex  functions  of  the  thermodynamic  quantities.  In  Section 
II  we  derive  the  SGS  scalar  variance  equation  and  show  that  it  includes  a  large  number  of 
subgrid  terms  the  modeling  of  which  is  uncertain.  Considering  this  uncertainty,  we  next 
examine  in  Section  III  the  conserved  scalar  extracted  from  the  DNS  database  and  explore 
whether  it  can  be  represented  in  LES  by  a  presumed  PDF  for  which  the  SGS  scalar  variance 
would  be  needed.  Further,  in  Section  IV,  we  assess  the  ability  of  two  models  to  accurately 
portray  the  SGS  scalar  variance.  The  databases  are  summarized  in  Section  V.  Model 
assessments  are  here  performed  on  these  databases  representing  transitional  states  obtained 
from  simulating  mixing  of  two  chemical  species  under  supercritical  pressure.  The  database 
consists  of  three  sets  of  species  and  was  created  in  the  context  of  a  temporal  mixing  layer. 
The  mixing  situation  has  been  here  selected  rather  than  a  chemically  reacting  case  because 
the  former  represents  a  conservative  choice  in  that  the  scalar  could  have  a  considerably  more 
complex  distribution  in  a  reactive  flow  due  to  coupling  among  thermodynamic  variables 
that  would  introduce  increased  flow  structure,  so  that  if  a  model  is  not  accurate  for  the 
mixing  situation,  it  will  certainly  be  even  less  accurate  for  a  reactive  flow.  On  the  other 
hand,  an  acceptable  model  describing  the  conserved  scalar  only  during  mixing  may  not  be 
acceptable  for  reactive  flows  but  it  represents  the  departing  point  for  constructing  such  a 
model.  Additionally,  one  commonality  among  the  evolution  of  the  flow  for  all  layers  was  the 
formation  of  High  Density- Gradient  Magnitude  (HDGM)  regions  populating  the  entire  flow, 
which  originate  from  the  combination  of  the  distortion  of  the  original  boundary  separating 
the  two  fluids  and  the  mixing  of  species  having  disparate  molar  masses.  These  HDGM 
regions  are  the  sites  where  mixing  between  species  occurs.  Thus,  the  conserved  scalar, 
which  is  any  of  the  two  species  in  a  layer,  displays  high  non-uniformities  and  the  HDGM 
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regions  resemble  in  geometry  the  flamelets  much  studied  in  reactive  flows;  this  resemblance 
adds  relevance  to  the  databases  utilized  in  this  study.  Results  are  presented  in  Section  VI. 
A  summary  and  conclusions  are  offered  in  Section  VII. 


II.  THE  SGS  SCALAR  VARIANCE  TRANSPORT  EQUATION 


Considering  Z  to  be  one  of  the  species,  the  SGS  scalar  variance  is 

az={Z2j-{Z}2  (1) 


where  the  operator  {*}  denotes  Favre-filt-ered  quantities  (we  depart  in  this  Section  from 
the  typical  Z  notation  given  the  complex  expressions  in  the  SGS  scalar  variance  equation 
below),  i.e. 

-pZ_ 

P 


where  p  is  the  mixture  density  and  for  any  function  g(x)  the  filtered  value  is  expressed  in 
physical  space  by 

^(x)  =  [ fi,(x/)G'(x  -  X'MX'  (3) 

Jn 

where  G  is  the  filter  function  associated  with  the  characteristic  filter  size  A  and  Q  is  the 
entire  spatial  domain.  According  to  the  scalar  equation  under  supercritical  conditions  [8-12], 
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the  transport  equation  of  the  SGS  scalar  variance  is 

=  V  ’  \P  iDaD}  Vffz]  +  2 p  {DaD}  V  {Z}  ■  V  {Z} 


Dt 


term  1 


term  2 


2 p  [{Dap}  {VZ  •  VZ}}  -  +  2  {Zj  V  •  (pr/) 

term  4 


term  3 


term  5 


2 p  [{DaDVZ  ■  VZ}  -  { DaD }  {VZ  ■  VZ}] 


term  6 


+  2V  •  [p  {DaDZ  VZ}  -  p  {DaD}  {Z  VZ}} 


term  7 


-  2  {Z}  V  •  Ip  { DaDVZ }  -  p  {Dan}  { VZ}] 


term  8 


2V  •  p  —  Z)Z2—  Vr|  —  2p  |«ba'(1  —  Z)Z— VT  ■  Vzj 


term  9 


term  10 


2V 


f  (1  —  Z)Z2  17111712  .  ^ 

p\D{—pr tp - Vp 


FtnT 


m 


term  11 


RVT 


m 

term  12 


—  2  {Z}  V  • 


D 


P  \  &bk(  1  —  Z)Z—  VT 


—  2  {Z}  V  • 


P{D 


term  13 

(1  —  Z)Z  mim2 


R„T 


m 


A  Vp 


(4) 


term  14 


where  ^  +  {u}  •  Vcrz  is  the  material  derivative,  t  is  time,  u  denotes  the  velocity, 

p  is  the  pressure,  T  is  the  temperature,  Ru  is  the  universal  gas  constant,  D  and  an  are 

the  variable  diffusion  coefficient  and  mass  diffusion  factor  respectively,  anx  is  the  Bearman- 

Kirkwood  form  of  the  thermal  diffusion  factor  and  A  is  defined  as 

/I  d(m/p)  1  d{m/p)\  (  ) 

\m2  dX2  mi  3Xl  )  ’  U 

with  rn  denoting  the  molar  mass  and  X  labeling  the  molar  fraction  (for  both  m  and  X, 

subscript  2  refers  to  Z  and  subscript  1  refers  to  the  complement  of  Z  in  the  mixture, 

(1  —  Z)).  Vectors  u  and  r/  are  the  SGS  fluxes 


u 


=  (Z2 u}  -  { Z2 }  {u}  , 


(6) 
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7]  =  {Zu}  -  {Z}  {u}  .  (7) 

Equation  (4)  clearly  differs  from  the  classical  SGS-scalar-variance  transport  equation 
under  atmospheric  conditions  (e.g.  Jimenez  et  al.  [6]).  On  the  right  hand  side  (r.h.s.) 
of  Eq.  (4),  one  recognizes  familiar  terms  such  as  the  molecular  diffusion  (term  1),  the 
Fick-issued  resolved  and  filtered  scalar  dissipation  rate  (terms  2  and  3,  respectively),  the 
transport  of  the  square  of  the  scalar  and  of  the  scalar  by  the  SGS  turbulence  (terms  4  and  5, 
respectively);  however,  new  terms  appear  as  a  consequence  of  the  spatial  variation  of  Dap 
under  supercritical  conditions  and  also  because  of  filtering  the  nonlinear  pressure-gradient 
and  Soret  terms.  We  distinguish  here  between  the  total  dissipation,  which  is  the  irreversible 
entropy  production  [13] 

X  oc  Ja  •  JQ,  (8) 

and  which  is  the  sum  of  six  terms  (since  the  species  mass-diffusion  flux,  JQ,,  is  the  sum  of 
three  terms  [8,  9,  13]),  and  the  Fick-issued  dissipation  which  represents  only  one  of  these 
six  terms. 

Whereas  the  form  of  Eq.  (4)  highlights  the  new  terms  due  to  the  supercritical-pressure  as¬ 
pect,  another  form  of  the  SGS  scalar  variance  equation  can  highlight  new  SGS  contributions 
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to  the  SGS  scalar  variance;  this  latter  form  is  similar  to  that  of  Pera  et  al.  [7]: 

P~r 7^=  V-lT>{DaD}Voz]  -V-[Zf(w-2{Z}i,)] 

Jut  s - V - '  s - V - ' 

Resolved  Molecular  Diffusion  SGS  Turbulent  Fluxes 

-  2 p  {Dap}  [{VZ  •  VZ}  ~V{Z}-V  {Z}\  -  2 pq  •  V  {Z} 

Fick-issued  SGS  Dissipation  SGS  Production 

+  2V  •  [; d{PaDZVZ }  -p{PaD}  {ZVZ}\ 

SGS  Diffusivity  1 

-  2' V  •  [p  {Z}  ({ DaDVZ }  -  {DaD}  V  {Z})\ 

" - v - ' 

SGS  Diffusivity  2 

+  2  [p  {DaDWZ}  -  p  {DaD}  V  {Z}]  •  V  {Z} 

" - V - ' 

SGS  Diffusivity  3 


2 p  [{DaDVZ  ■  VZ}  -  {DaDj  {VZ  ■  VZ}] 


2  V  • 


SGS  Diffusivity  f 

,D. 
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SGS  Soret  1 


2  p 


aBK(l  -  Z)Z^VT  ■  Vz\  -  |  aB*r(l  -  Z)Z^VT  }-V{Z} 


SGS  Soret  2 


■  2  V 


/  ,  „  (1  —  Z)Z2  777-17712  .  _ 

P  D—rrP: - —A  Vp 


RuT 


m 


c(i \  {z } 

RV1  m 


SGS  Pressure  1 


-2  p 


D(l-Z)Zmm,AVp .  vZ)  _  ( A  .  v  {z} 


RVT  m 


R„T  m 


(9) 


SGS  Pressure  2 


where  term  5  in  Eq.  (4)  was  split  into  the  classical  turbulent-flux  form  and  SGS  production, 
the  turbulent  fluxes  were  combined  and  only  SGS-t-ype  terms  emphasized.  The  meaning  of 
the  new  SGS  terms  is  not  obvious  and  must  be  investigated. 

It  will  be  shown  in  Section  VI A  that  the  additional  SGS  terms  in  Eq.  (4)  have  significant- 
activity,  and  thus  are  important  to  model.  Moreover,  the  root  mean  square  (r.m.s.)  activity 
and  the  mean  values  of  the  SGS  terms  in  Eq.  (9)  will  be  examined  to  understand  their 
influence  on  the  SGS  variance  evolution.  Models  for  the  new  SGS  terms  in  either  Eq.  (4) 
or  Eq.  (9)  are  not  currently  available,  a  fact  which  motivated  investigating  whether  the 
conserved  scalar  conforms  to  a  PDF  having  a  known  mathematical  expression.  If  this  is  the 
case,  the  PDF  could  be  used  in  a  presumed-PDF  approach  and  therefore  the  SGS  scalar 
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variance  could  be  useful  for  constructing  the  PDF,  a  fact  which  would  motivate  finding  a 
model  for  it.  The  investigation  of  the  SGS  PDF  shape  is  described  next.  In  the  remaining 
part  of  this  study,  we  revert  to  the  usual  notation,  Z,  for  the  Favre  filtered  quantities. 


III.  THE  PDF  OF  THE  CONSERVED  SCALAR 

The  scalar  statistics  is  the  important  information  desired  for  flamelet  modeling  (e.g. 
[14]).  However,  if  the  PDF  is  complex,  it  has  the  drawback  that  it  can  only  be  approximately 
reconstructed  by  a  large  number  of  its  moments  [15].  It  is  certainly  computationally  easier  if 
it  can  be  shown  that  the  PDF  of  the  scalar  conforms  to  a  PDF  having  a  simple  mathematical 
form  for  which  only  a  small  number  of  moments  is  necessary  for  its  reconstruction.  Such 
simple  PDFs  are,  for  example,  the  Dirac,  the  Gaussian  and  the  /?  functions  for  which  one 
needs  a  maximum  number  of  two  moments  for  reconstruction.  As  stated  above,  the  first 
moment  of  the  local  (SGS)  PDF,  is  always  computable  from  LES,  so  that  the  focus  is  on 
computing  the  second  moment,  i.e.  the  SGS  scalar  variance.  An  extensive  literature  exists  on 
the  topic  of  the  SGS  scalar  variance  computation  [2,  4,  5,  7,  16,  17]  but  not  for  supercritical 
turbulent  flows  for  which,  as  seen  in  Section  II,  the  variance  equation  is  considerably  more 
complex  than  in  the  already  studied  flows,  indicating  possible  departures  from  previous 
findings,  as  implied  by  previous  results  [8].  To  inquire  about  the  form  of  the  SGS  PDF, 
we  recall  here  its  definition  and  discuss  a  practical  way  to  calculate  it  under  compressible 
conditions. 

The  instantaneous  SGS  PDF  (fsgs)  of  a  scalar  quantity  Z  may  be  defined  using  the  fine¬ 
grained  density  y[£,  Z(x)]  =  c>[Z(x)  —  £],  by  weighting  it  with  the  filter  function  G  [18,  19] 
as 

fsgsit;  x)  =  f  £[Z(x')  -  £]G(x  -  x)c/x  (10) 

Jn 

where  5  is  the  Dirac  delta  function  and  £  is  the  scalar-space  variable  representing  Z.  If  G  is  a 
positive  filter  kernel,  fsgs  has  the  property  of  a  PDF.  As  every  PDF,  it  is  constrained  by  the 
normalization  condition  stating  that  the  integral  over  the  scalar  space  is  unity.  For  Z  6  [0,1], 
the  instantaneous  filtered  value  of  any  quantity  g  (Z(x))  may  be  obtained  by  integration  over 
the  scalar  space  as  [17] 

g(x)  =  [  g(0fsgs(^x)d^.  (11) 

j  0 
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In  the  same  spirit,  for  obtaining  the  Favre- filtered  value  of  every  quantity  g  (Z(x))  by  using 
a  PDF,  a  mass-weighted  SGS  PDF  must  be  defined  [20] 

fsgsM  x)  =  =^y  ^p(x')^(x')  ~  £]G(x  -  x')dx'  (12) 

where  the  subscript  c  denotes  the  compressible  situation,  and  then 

ff(x)=  [  9(0fsgsM^)dC  (13) 

do 

The  mass-weighted  SGS  PDF  (referred  to  as  Filtered  Mass  Density  Function  (FMDF)) 
is  defined  by  Jaberi  et  al.  [21]  and  Raman  et  al.  [22]  in  a  probabilistic  manner  using  the 
property  that  the  integral  of  fsgSc  over  the  sample  space  gives  the  filtered  density  value.  This 
definition  is  then  used  in  the  framework  of  a  PDF  approach  where  a  transport  equation  for 
the  PDF  is  solved.  Here,  instead,  a  PDF  normalized  by  p(x),  i.e.  Eq.  (12),  is  needed  as  we 
wish  to  satisfy  the  classical  definition  of  the  PDF  integrating  to  unity.  The  SGS  PDF,  as 
defined  by  Eqs.  (10)  and  (12),  is  not  a  statistical  quantity  as  it  is  a  one-point  distribution 
conditioned  to  a  given  flow  realization.  Following  Jimenez  et  al.  [3],  Fox  [23]  and  Pitsch  [24], 
an  appropriate  SGS  PDF  utilizable  for  modeling  purposes  should  instead  be  understood  as 
a  statistical  quantity  arising  from  a  statistical  sample  of  equivalent  grid  elements.  Since  in 
the  a  priori  analysis  using  the  DNS  database  the  exact  statistical  quantity  is  not  available, 
we  will  instead  use,  as  an  approximation  of  the  statistical  SGS  PDF  for  evaluating  the 
accuracy  of  models,  statistics  computed  employing  the  one-point,  SGS  PDF  defined  above. 
The  expectation  of  this  PDF  over  homogeneous  directions  and  the  conditional  (on  moments) 
expectation  are  used  to  assess  presumed  PDF  shapes  against  the  DNS-extracted  ones,  as 
described  below. 

Assuming  the  filter  width  to  be  smaller  than  the  variation  scale  of  mean  quantities  and 
that,  filtering  and  averaging  operators  commute,  a  Favre  mean  (i.e.  expectation)  of  the 
filtered  local  value  may  be  obtained 

<  <7(x)  >c=  [  g(0  <  fsgsc(&x)  >c  d£  (14) 

Jo 

where  the  operator  <  •  >c  denotes  planar  averages  weighted  by  the  filtered  density.  The 
quantity  <  fsgSc  >c  is  equivalent  to  a  filtered  PDF  of  the  scalar  as  defined  by  Jimenez  et 
al.  [3]  for  incompressible  flows  and  it  is  given  by  the  mass-weighted  PDF  of  the  scalar  over 
planes  of  height  equal  t-o  the  filter  width  (coarsened-grid  planes).  The  filtered  PDF  may  be 
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computed  as 


<  fagac{&  G  >c=  J  /c(G  V^fagadb  G  0*)dfefo*  (15) 

by  using  the  mass- weighted  joint-moment  PDF  fc  in  conjunction  with  the  PDF  fsgSc,  defined 
as  a  function  of  couples  of  moments  (G  ay)  instead  than  x  where  ay  is  the  moment-space 
variable  representing  aZ-  and  integrating  over  the  LES-moment  space  (G  ay).  The  defini¬ 
tion  of  Eq.  (15)  may  be  computationally  practical  since  any  assumed  PDF  model  may  be 
evaluated  by  replacing  in  Eq.  (15)  the  presumed  PDF,  /S9Sc(G  £,  cr^),  computed  from  any 
couple  of  mapped  values  (G  ay),  and  /C(G  oy)  as  obtained  from  the  filtered  DNS  database 
over  homogeneous  (coarsened-grid)  directions.  The  fsgSc  expectation,  <  fsgSc  >c,  is  used  to 
evaluate  the  presumed-PDF  shapes  against  the  DNS-extracted  ones. 

To  approximate  the  statistical  SGS  PDF,  we  also  consider  as  an  alternate  to  the  above 
method,  the  optimal  estimator  method  of  Moreau  et  al.  [25]  which  uses  the  conditional  ex¬ 
pectation.  The  latter  is  obtained  by  averaging  the  PDF  fsgSc(£',  G  ay),  over  a  sample  of  PDFs 
having  the  same  couples  of  exact  moments  (G  ay).  Practically,  the  computation  is  performed 
by  averaging  over  SGS  volumes  having  couple  of  moments  (G  a^)  very  close  to  a  selected  one, 
according  with  a  fixed  standard  deviation  chosen  such  as  to  have  accurate  statistics.  The 
conditional  expectation  of  the  DNS-extracted  quantity  /sSSc(G  G  &$)  is  then  compared  to  the 
assumed  PDFs  computed  using  the  exact  moments.  Here,  the  Favre-conditional  expecta¬ 
tion  is  obtained  by  mass-weighted  averaging  using  the  filtered  local  density,  and  statistics 
are  computed  over  a  slab  of  the  mixing  layer  in  order  to  increase  the  sample  size. 

In  this  study,  three  local  PDFs  will  be  assessed:  the  (3  PDF,  the  Gaussian  and  the  Dirac 
ones.  These  PDFs  are  here  briefly  recalled: 


The  (3  PDF: 


fsgSc  (g«(g^),g(g^))  =  1 


r(«)r(/i) 

where  F  is  the  Gamma  function  and  parameters  a  and  f3  are  defined  by 


a  =  {,«lrO_i 


The  Gaussian  PDF: 

fsgsc{£,  G  T;)  — 


0  =  (1-0  I  — — — -i 


J 


1  e  2as 


2l T<7c 


U1  +  erfoil 


i(1  +  “f:S: 


(16) 


(17) 


(18) 
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As  the  scalar  £  is  bounded  between  0  and  1,  the  Gaussian  distribution  is  truncated  and 
the  resulting  PDF  is  thus  re-normalized  using  the  difference  between  the  cumulative 
distribution  functions  evaluated  at  b  =  1  and  a  =  0.  This  computed  distribution  has 
no  longer  the  same  mean  and  variance  as  the  original  distribution.  As  an  alternate,  a 
clipped  Gaussian  is  also  considered  [26]. 

•  The  Dirac  PDF: 


f, „.((;()  =  <%-«•  (is) 

The  Dirac  PDF  is  utilized  when  the  scalar  £  is  only  modeled  through  its  mean,  £. 

In  Section  VI B  we  present  results  from  computations  with  these  PDFs  using  both  meth¬ 
ods  to  approximate  the  statistical  SGS  PDF,  and  particularly  we  show  that  whereas  the 
Dirac  and  Gaussian  PDF  are  generally  deficient,  the  /5-PDF  typically  yields  a  good  ap¬ 
proximation,  a  fact  which  provided  the  incentive  to  directly  model  the  scalar  variance,  as 
described  in  the  following. 


IV.  DIRECT  MODELING  OF  THE  SGS  SCALAR  VARIANCE 

Two  types  of  models  appear  to  be  promising  candidates  for  modeling  small-scale  effects 
removed  through  filtering.  The  first  model  is  of  a  structural  type  [27]  based  on  the  Ap¬ 
proximate  Deconvolution  Model  (ADM)  [28]  while  the  second  one  is  of  a  functional  type 
using  a  gradient-based  scaling  law  [4] .  The  two  models  are  conceptually  different  as  the  first 
one  arises  from  a  mathematical  derivation  with  no  assumption  regarding  the  nature  of  scale 
interactions,  while  the  second  one  uses  a  mixing- length  hypothesis  in  conjunction  with  an 
equilibrium  assumption.  The  original  contribution  of  our  work  is  the  investigation  of  the 
ADM  capability  for  computing  az  in  variable  density  flows,  and  the  extension  of  a  recent 
dynamic  gradient-based  formulation  [5]  to  compressible  conditions  for  modeling  az- 
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A.  The  approximate  deconvolution  model 


The  deconvolution  procedure  [28]  relies  on  the  assumption  that  there  exists  an  inverse 
operator  GG1  =  YliL o(I  ~  G)1  such  as 

where  I  is  the  identity  operator  and  Z*  represents  an  approximation  of  the  original  field  Z, 
on  the  LES  mesh  grid,  at  the  series  truncation  order  N.  The  series  is  known  to  converge  for 
|| /  —  6*11  <  1,  and  written  using  de  van  Cittert  method  [29]  based  on  the  Neumann  series 
leads  to 

Z*  =  Z+  (z-f)  +  +  + .  (21) 

where  the  accuracy  of  Z*  depends  on  N .  In  practice,  due  to  the  numerical  discretization 
(projection  on  a  coarse  mesh  grid),  only  the  recoverable  part  of  the  original  field  Z  may  be 
obtained  by  a  deconvolution  procedure  [30].  Furthermore,  the  deconvolution  relies  on  the 
form  of  the  convolution-filter  kernel  and  on  the  shape  of  the  spectrum  of  the  field  [17]. 

The  series’  rate  of  convergence  is  situation  dependent.  For  example,  it  has  been  found 
[31]  that  N  =  3  is  sufficient  to  bring  an  improvement  in  that  for  a  quantity  q i>,  ((f)*  —  q !>)  is  not 
null  showing  that  the  expansion  is  indeed  effective,  and  for  N  >  5  the  value  of  ((f)*  —  o)  did 
not  change  appreciably  from  that  obtained  with  N  <  5.  In  other  studies  specifically  directed 
at  Z,  it  was  found  that  the  series  converges  very  slowly.  Pantano  and  Sarkar  [16]  tested 
the  deconvolution  procedure  in  an  a  priori  analysis  using  the  DNS  of  a  temporal  mixing 
layer;  they  showed  that  even  with  N  =  5  (i.e.  fourth-order  approximation)  and  a  small 
filter  size  (A/A x^ns  =  4,  where  A xdns  is  the  grid  spacing  for  a  simulation  where  all  scales 
relevant  to  most  of  the  dissipation  are  resolved,  as  in  DNS),  no  more  than  88%  of  the  total 
SGS-scalar-variance  amount  (in  the  peak  zone)  was  recovered.  The  rate  of  convergence  of 
the  series  and  the  recoverable  amount  of  the  field  are  thus  the  main  issues  of  the  method. 
In  order  to  recover  all  effects  of  the  smallest  scales,  Stolz  et  al.  [28]  suggested  a  secondary 
filtering  through  the  use  of  a  relaxation  parameter  in  the  Navier-Stokes  (NS)  equations.  To 
mitigate  both  these  issues,  Pantano  and  Sarkar  [16]  and  Mellado  et  al.  [17]  proposed  an 
approximate  reconstruction  using  moments  (ARM)  method  which,  based  on  the  definition 
of  an  intermediate  field,  leads  to  computing  the  SGS  scalar  variance  by  imposing  equality 
between  the  real  moment  and  that  of  the  intermediate  field.  This  procedure  results  in  an 


£(/  -  G)‘  *z  <2°) 
1=0 
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expression  for  a  parameter  Co,  characterizing  the  model,  that  is  obtained  by  finding  the  real 
and  positive  solution  of  a  second-order  polynomial  equation.  The  polynomial  coefficients 
a0,  ai,  a 2  are  computed  using  the  presumed  shape  of  the  spectrum  of  the  scalar  quantity. 
Once  the  shape  of  the  spectrum  is  selected,  c0  may  be  precomputed  and  stored  in  a  two- 
dimensional  table  of  coordinates  represented  by  two  dimensionless  quantities:  the  Peclet 
number,  and  the  ratio  between  the  filter  width  and  a  large  scale  of  the  scalar  fluctuations. 
That  approach  was  developed  for  incompressible  flows.  Under  compressible  conditions,  that 
model  requires  rewriting  for  the  conservative  quantities  for  which  the  governing  equations 
are  solved  in  LES  of  compressible  flows;  but  then,  equality  between  moments  would  not 
lead  to  an  explicit  expression  for  the  SGS  scalar  variance  unless  additional  assumptions  are 
invoked.  Moreover,  when  in  presence  of  variable  density,  the  shape  of  the  spectrum  cannot 
be  built  using  calibrated  constants  computed  either  by  using  the  Reynolds-averaged  scalar 
variance  or  by  employing  the  Favre-averaged  one.  This  is  why,  considering  the  practical 
utilization  in  LES,  we  opt  to  explore  a  classical  deconvolution  procedure. 

Since  in  LES  it  is  the  conservative  rather  than  primitive  quantities  which  are  calculated, 
we  are  legitimately  inquiring  as  to  whether  the  SGS  scalar  variance  could  be  accurately 
reconstructed  employing  a  deconvolution  procedure  for  these  quantities.  Indeed,  for  recon¬ 
structing  the  approximated  scalar  field  and  building  the  SGS  (Favre)  scalar  variance,  we 


must  reconstruct  p  and  Z  through 

P*  =  P  +  {p  -  p)  +  (p  -  2P  +  d)  + .  (22) 

(pZ)  =  pZ  +  ^ pZ  —  pZ^j  +  ^ pZ  —  2  pZ  +  pZ^j  + .  (23) 

which  then  permits  writing 

Z**  =  (pZ)* /p*.  (24) 

invoking  the  assumption 

Z**  ~  Z*.  (25) 


Equation  (25)  implies  equality  between  second-order  moments  of  Z**  and  Z* .  The  non¬ 
linear  function  of  the  approximate  field,  the  Favre  SGS  scalar  variance,  is  then  built,  for 
consistency  [16],  by  using  of  deconvoluted  fields  as 

Z  p*  p*  p* 
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The  assumption  of  Eq.  (25)  is  exact  if  both  p  and  pZ  fields  are  totally  recovered  and  an 
infinite  expansion  is  used.  Otherwise,  there  is  no  proof  that  ratio  between  series  (23)  and 
(22)  converges  to  the  recoverable  part  of  the  quantity  Z .  which  is  Z*  at  a  given  truncation 
order  N.  Such  an  assumption,  although  not  trivial,  has  already  been  used  by  Stolz  et  al. 

[32]  while  reconstructing  flux  terms  of  the  compressible  NS  equations,  and  by  Dubois  et  al. 

[33]  in  the  first  of  a  two-step  procedure  for  estimating  the  SGS  stress  tensor.  The  difference 
between  the  previous  studies  and  our  investigation  is  that  ADM  is  here  used  in  the  context  of 
a  soft-deconvolution  problem  ([27])  in  which  no  additional  models  for  recovering  the  deficient 
SGS  part  are  provided.  In  Section  VI C  1  we  evaluate  the  assumption  of  Eq.  (25)  and  the 
equality  between  moments  of  Z**  and  Z*  for  several  truncation  orders.  We  also  assess  the 
ADM  model  using  the  series  Eqs.  (22)  and  (23),  with,  consistently,  both  series  truncated  at 
same  order. 

B.  The  dynamic  model 

A  Smagorinsky-type  model  for  predicting  the  SGS  scalar  variance  of  compressible  flows 
was  proposed  by  Pierce  and  Moin  [4]  using  a  dynamic  procedure  [34],  This  model  was 
similar  to  that  of  Moin  et  al.  [35]  for  modeling  the  SGS  stress  tensor  and  heat  flux  under 
compressible  conditions.  The  model  of  Pierce  and  Moin  [4]  has  been  extensively  used  over  the 
years  [7,  20,  36,  37]  as  an  alternate  to  the  scale-similarity  model  suggested  by  Cook  and  Riley 
[2]  based  on  the  idea  of  Bardina  et  al.  [38] .  The  scale-similarity  model  has  been  explored  in 
several  different  configurations  [3,  7,  39,  40].  Recently,  Balarac  et  al.  [5],  using  the  optimal 
estimator  concept  of  Moreau  et  al.  [25] ,  showed  that  the  irreducible  error  associated  with 
a  Smagorinsky-type  model  is  relatively  small  if  compared  to  that  evaluated  for  a  scale- 
similarity  model,  meaning  that  the  functional  form  of  a  gradient-based  model  has  significant 
potential  in  LES  for  reproducing  the  SGS  scalar  variance.  This  has  also  been  observed  by 
Wall  et  al.  [36]  when  studying  the  performance  of  the  two  models  used  in  that  study  in 
conjunction  with  a  presumed  PDF  approach.  On  the  other  hand,  Balarac  et  al.  [5]  showed 
that  the  quadratic  error  associated  with  the  gradient-based  model  notably  increases  with 
filter  width  when  a  classical  dynamic  procedure  is  used  for  computing  the  model  coefficient. 
In  the  present  study,  the  Balarac  et  al.  [5]  model,  which  is  for  incompressible  flow,  is 
reformulated  under  compressible  conditions  and  evaluated  for  compressible  supercritical 
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turbulent  mixing  layers  in  Section  VIC 2. 

For  completeness,  the  Pierce  and  Moin  [4]  dynamic  model  is  briefly  recalled  and  for 
consistency  we  use  their  notation,  where  Q  is  the  unweighted  test  filter  at  test-filter  level  G 
corresponding  to  the  filter  width  A  and 

z  =  izfp-  (27) 


In  the  classical  dynamic  model  [4] ,  first  the  variance  is  related  to  the  gradient  of  the  resolved 
scalar  field  as 


~poz  =  p(ZZ  -  ZZ)  =  Cd A  p 

which,  when  filtered,  yields  the  following  expression 


vz 


(28) 


pZZ  -  pZZ  =  Cd A  p 


VZ 


(29) 


A  similar  model  is  then  proposed  at  the  filter  level  corresponding  to  the  double  convolution 
of  G  and  G,  i.e.  G  =  G  *  G,  associated  to  the  filter  width  A  : 

. "'"V  —  2  i  8^  1 2 


pZZ  -  pZZ  =  CdA  p 


VZ 


(30) 


Finally,  subtracting  Eq.  (29)  from  Eq.  (30)  and  assuming  C'd  =  G”  =  Cd  (implying  either 
that-  the  flow  variation  is  such  that  filtering  with  the  grid  or  test  filter  produces  analogous 
fields  or  that  the  coefficients  have  a  slow  variation)  yields 


_ 

^2^ 

2  2  / 

~~ 

2\ 

pZZ  -  pZZ  =  cd 

A  p 

VZ 

-A2(p 

VZ 

) 

(31) 


which  may  be  generically  written  as 


C  =  CdM 


where  C  = 


\ 


■pZZ  -  -pZZ 


V 


/ 


is  the  generalized  Leonard  term  and  Ai  is  defined  as 


(32) 


M 


^2^ 
A  p 


VZ 


(33) 


Since  both  £  and  Ai  are  computable  from  the  LES  solution,  Cd  may  be  dynamically  com¬ 
puted  over  homogeneous  directions  either  as  the  ratio  between  averaged  £  and  Ai,  or  through 
the  least-square  averaging  (Lilly  [41])  which  optimizes  the  local  value  by  minimizing  the 
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quadratic  error.  This  dynamic  procedure  uses  a  compressible  version  of  the  Germano  iden¬ 
tity  [42],  originally  written  for  the  SGS  stress  tensor  [35,  43] 

h  ij  P^ij  P^~ij  >  (34) 

which  adapted  to  the  scalar  variance  may  be  written  as 

L  =  pEz  -  pcrz  (35) 


where 

%z  =  (pZZ  -  .  (36) 

Equation  (30)  was  questioned  by  Balarac  et  al.  [5]  who  emphasized  that  a  generic  fil¬ 
tered  quantity  /  is  not  directly  obtained  through  a  single  convolution,  but  rather  from  two 
sequential  filterings.  Using  the  Bedford  and  Yeo  [44]  expansion  (which  provides  a  power 
series  for  the  non-linear  filtered  generic  term  fg  as  a  function  of  the  resolved  quantities  / 
and  g ),  applied  twice  to  the  filtered  quantities,  Balarac  et  al.  [5]  pointed  out  a  missing 
leading-order  term  in  Eq.  (30)  and  provided  an  alternative  model  for  Ai  that  involves  the 
use  of  the  Leonard  term  expansion.  The  same  reasoning  is  possible  for  compressible  flows 
using  an  appropriate  Taylor  expansion.  Such  an  approximate  expansion  was  suggested  by 
Vreman  [43]  for  modeling  the  SGS  stress  tensor  and  is  adapted  here  for  isotropic  filter  to 
the  scalar  variance  as 


pG z  =  pZZ  —  pZ  pZ  j p  (37) 

X2  X2  X2  X2 

+  -V2  (pZZ) -{pZ  +  -V2  (PZ))(pZ  +  -V2  (pZ))/(p  +  — V2p)  +  0(A4) 
X2 

=—p\vz\2  +  o(A4) 


where  the  approximation 


A2  , 


1  A  V2p  +  0(A4) 


p  24  p2 

has  been  made  to  obtain  the  final  result.  A  similar  expansion 


(38) 


p  =  p+0(  A2),  Z  =  Z+0(A2) 


(39) 


used  in  Eq.  (37)  leads  to 


pZZ 


_ —  A  _ 

pZZ  +  —p 


vz 


-r-4 , 


0(  A’). 


(40) 
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Equation  (40)  is  conceptually  similar  to  that  used  by  Balarac  et  al.  [5]  for  incompressible 
conditions  and  in  the  same  way  it  may  be  used  for  obtaining  i.e.  the  left  hand  side 
(l.h.s.)  term  of  Eq.  (30).  First,  the  test  filter  G,  is  applied  to  Eq.  (40) 


pZZ  =  pZZ  +  —  (  P 


vz 


0(  A4), 


(41) 


then  the  first  term  on  the  r.h.s.  is  expanded  at  the  filter  level  G  based  on  p  and  Z  using  the 
approximate  expansion  of  Eqs.  (37)-(39)  applied  to  these  variables,  leading  to: 


— -~~y  ^ 

2  A  / 

_ _ _ 

pZZ  =  J)ZZ  +  —  p 

VZ 

+  w 

VZ 

+  0(A  ,  A  ). 


(42) 


The  above  equation,  reformulated  using  a  constant  to  account  for  the  truncation  error,  leads 
to 


pZZ  -  pZZ  =  Cd 


A2^ 


vz 


A2 


P 


vz 


(43) 


which  represents  another  form  for  the  l.h.s.  term  in  Eq.  (30).  Subtracting  Eq.  (29)  from 
Eq.  (43)  and  assuming  C'd  =  Cd  =  Cd  (under  the  same  implications  as  when  adopting 
C'd  =  =  Cd)  leads  to  a  new  formulation  of  At,  Mn 

Mn  =  A2J>  VZ  2  (44) 


representing  the  leading  order  of  the  Leonard  term  expansion  when  C,i  palliates  for  the 
truncated  terms  in  the  Taylor  series  expansion.  Balarac  et  al.  [5]  pointed  out  that  Eq.  (42) 
contains  a  new  leading  order  term  with  respect  to  Eq.  (30)  which,  when  taken  into  account, 
yields  the  new  formulation  of  At,  Atn.  Although  we  obtain  the  same  result  as  Balarac 
et  al.  [5],  Mn  for  compressible  conditions,  Eq.  (30)  is  different  from  the  corresponding 
equation  of  Balarac  et  al.  [5]  as  far  as  filter  width  used.  Balarac  et  al.  [5]  did  not  discuss 
the  disparity  between  the  filter  width  used  in  the  classical  formulation  and  the  filter  width 
needed  in  the  Leonard  term  expansion,  and  it  appears  that  the  same  filter  width  was  used. 
However,  when  the  correct  filter  width  is  used,  the  issue  of  the  new  leading  order  term  in 
the  formulation  of  Eq.  (30)  is  moot,  and  the  result  is  that  only  an  alternative  formulation 
for  J>Zz  is  derived.  According  to  Eq.  (44),  it  is  clear  that  the  new  formulation  uses  A.  For 
the  classical  formulation,  according  to  Eq.  (33),  the  filter  is  though  A.  This  finding  is  in 
agreement  with  the  study  of  Brun  and  Friedrich  [45] .  Using  the  Vreman  et  al.  [46]  estimate, 
the  filter  width  associated  with  the  filter  level  G  is 

A2  =  A2+A2,  (45) 
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representing  an  approximation  for  double  top-hat  filtering.  In  order  to  explore  the  impact 
of  the  filter  width  used  in  the  classical  formulation,  and  thus  to  asses  the  correctness  of  the 
results,  we  present  in  Section  VI C  2  an  a  priori  analysis  using  the  classical  formulation  in 
conjunction  with  each  of  the  approximations  for  the  filter  width  A  at  the  filter  level  G:  i) 
A  =  A2  +  A2  (correct  approximation)  and  ii)  A  =  A  (incorrect  with  a  top-hat-  filter), 
and  compare  the  results  from  these  two  computations  with  those  from  the  new  formulation 
based  on  the  Leonard  term  expansion  for  which  the  filter  is  demonstrated  to  be  A. 

We  can  thus  evaluate  the  gradient-based  model  for  az  using  two  different  formulations. 
The  first  formulation  is  the  classical  one  employing 

(CM) 

d  (MM) 

where  ()  denote  averages  over  homogeneous  planes.  The  second  formulation,  which  we  call 
the  “new  model”,  is  the  present  one,  and  uses 


{CMn) 

{. M„Mn >' 


(47) 


Since  the  evaluation  of  the  correct  filter  width  at  filter-level  G  represents  the  main  issue  for 
computing  Gcj  using  a  dynamic  procedure,  the  model  based  on  the  Leonard  term  expansion 
represents  a  solution  to  this  quandary,  as  shown  in  Section  VI C  2. 


V.  DESCRIPTION  OF  THE  DNS  DATABASE 

A  complete  and  detailed  description  of  the  DNS  database  has  already  been  provided  in 
[8-12],  Out  of  the  complete  database  consisting  of  three  sets  of  species  heptane/nitrogen 
(HN),  oxygen/hydrogen  (OH)  and  oxygen/helium  (OHe)  for  which  several  simulations  were 
conducted,  we  select  here  three  simulations,  HN600,  OH750  and  OHeGOO  for  examination 
so  as  to  enhance  the  generality  of  the  results.  The  DNS  were  conducted  for  a  temporal 
mixing  layer  and  initiated  with  (vorticity-thickness  based)  Reynolds  number  of  600  (HN600, 
OHeGOO)  and  750  (OH750)  where  8Ujq  =  A U0/  (dui/dx2)  is  the  initial  vorticity  thickness, 
()  is  here  performed  over  (x\.  x:>)  planes  and  A(/0  is  the  initial  velocity  difference  across  the 
layer.  In  all  cases,  the  DNS  grid,  Axdns,  was  fine  enough  to  resolve  the  scales  overwhelm¬ 
ingly  responsible  for  the  dissipation.  Transitional  states  were  achieved  for  all  of  these  layers 
at  t\r  =  135  for  HN600,  t*r  =  150  for  OH750,  t*r  =  220  for  OHeGOO,  where  t*  =  tAU0/5 W)0; 
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the  corresponding  transitional  Reynolds  numbers  were  1452,  1507  and  2004.  Selle  et  al.  [9] 
showed  that  this  database  is  relevant  for  fully-turbulent  flow  modeling. 


VI.  RESULTS 

In  all  computations,  the  results  of  which  are  presented  below,  we  use  a  top-hat  filter  and 
a  trapezoidal  integration  method. 

A.  Activity  of  terms  in  the  scalar  variance  equation 

All  terms  in  Eqs.  (4)  and  (9)  were  a  priori  evaluated  using  the  three  DNS  realizations 
discussed  above.  The  goal  of  this  evaluation  is:  i)  to  assess  the  magnitude  of  the  filtered 
terms  with  respect  to  that  of  the  resolved  terms  for  modeling  purposes  (Eq.  4)  and  ii)  to 
assess  the  magnitude  of  the  SGS  terms  with  respect  to  that  of  the  resolved  terms  in  order 
to  identify  which  SGS  quantities  are  important  to  model  (Eq.  9).  A  cubic  top-hat  filter  and 
several  filter  widths,  A,  were  used  in  the  evaluation  but  only  the  analysis  corresponding  to 
the  ratio  A/ Ax dns  =  8  is  here  presented,  being  representative  of  all  ratio  values. 

Depicted  in  the  Figure  (1)  is  the  activity  of  terms  in  Eq.  (4)  as  measured  by  the  r.m.s. 
magnitude  of  each  term.  The  results  indicate  that  terms  4  and  5,  which  involve  the  classical 
subgrid  fluxes,  dominate  showing  that  our  transitional  databases  have  indeed  the  turbulent 
characteristics  that  make  them  relevant  to  this  study.  The  importance  of  the  advective  term 
increases  with  the  strength  of  the  HDGM  at  transition  [11]  (compare  OH750  with  HN600) 
and  with  increasing  Reynolds  number  value  at  transition  [9]  (compare  HN600  and  OHe600). 
The  diffusion  term  (term  1)  is  of  same  order  of  magnitude  as  the  advect-ion  term  only  for 
OH750  because  the  hydrogen  diffusivity  is  very  large.  Terms  2  and  3  are  somewhat  smaller 
than  the  advection  term  for  HN600  and  OHeGOO,  but  significantly  larger  for  OH750.  For 
both  HN600  and  OHeGOO,  and  less  so  for  OH750,  term  3  is  larger  than  term  2,  which  is  a 
manifestation  of  the  SGS  magnitude  associated  with  the  Fick-issued  dissipation.  Term  6, 
representing  SGS  effects,  could  be  neglected  in  all  cases.  Evaluation  of  term  7,  expressing  the 
diffusivity  spatial  variation,  shows  that  it  is  of  the  same  order  of  magnitude  as  the  resolved 
diffusion  (term  1)  for  HN600  and  OHeGOO  simulations,  and  even  larger  for  the  OH750.  Thus, 
term  7  is  non  negligible  in  Eq.  (4)  under  all  circumstances.  Similarly,  term  8  representing 
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the  effect  of  spatially  varying  diffusivity  is  comparable  to  term  1  for  all  cases  and  should  be 
included  when  modeling  Eq.  (4),  particularly  when  diffusion  plays  an  important  role  during 
mixing  (i.e.  OH750).  In  contrast,  the  terms  arising  from  filtering  of  the  Soret  contribution, 
i.e.  terms  9,  10  and  13  are  only  substantial  for  the  HN600  case,  while  those  stemming  from 
filtering  of  the  term  proportional  to  the  pressure  gradient,  namely  terms  11  and  14,  only 
play  a  minor  role  in  the  OH750  and  OHeGOO  layers  and  are  unimportant  for  the  HN600  case. 
Term  12  is  found  negligible  everywhere. 

Figure  (2)  shows  the  results  from  the  a  priori  assessment  of  the  terms  in  Eq.  (9).  Both 
the  planar  average  and  the  r.m.s.  activity  are  computed  in  order  to  evaluate  not  only  the 
importance  but  also  the  nature  (i.e.  diffusive,  dissipative,  etc.)  of  each  term.  Considering 
the  planar  averages,  in  all  cases,  the  magnitude  of  the  diffusion  terms  is  comparable  to  that 
of  the  advection  term,  and  the  production  term  is  comparable  in  magnitude  to  the  Fick- 
issued  SGS  dissipation  term.  Due  to  the  difficulty  of  entraining  the  lower-stream  heavy  fluid, 
the  mixing  layer  growth  is  moderate  which  explains  the  small  value  of  the  advection  term 
compared  to  that  of  the  SGS  production  and  SGS  Fick-issued  dissipation.  Additionally,  for 
the  HN600  layer  the  term  denoted  by  SGS  Soret  2  has  considerable  negative  magnitude, 
adding  to  the  Fick-issued  dissipation.  This  observation  is  consistent  with  the  definition  of 
the  total  scalar  dissipation  [13] 

Xt  =  —r. - Jq  ■  J 

pDctD 

that  is  here  contrasted  to  the  Fick-issued  dissipation 

Xf  =  pDaoVZ-VZ,  (49) 

which  is  only  one  of  the  six  terms  of  Eq.  (48).  Obviously,  the  Soret  contribution  plays 
an  important  role  for  the  HN600  layer  in  the  destruction  of  the  scalar  fluctuations  at  the 
smallest  scales,  but  clearly  its  importance  is  binary-species  dependent  as  it  has  no  impact 
on  the  mixing  of  the  oxygen  with  hydrogen  or  helium.  In  order  to  understand  the  behavior 
of  xf  with  respect  to  xt ,  we  illustrate  in  Fig.  (3)  the  average  of  xr  conditioned  on  xf 
in  two  planes  of  the  HN600  mixing  layer  (a^/^.o  =  0.44  which  is  in  the  central  part  of 
the  layer  and  xg/fG.o  =  5.11  which  away  from  the  center  but  still  in  a  significant  mixing 
region).  The  results  indicate  that  xt  is  larger  than  xf,  particularly  at  the  periphery  of  the 
layer.  Noteworthy,  the  linear  dependency  between  xt  and  xf  in  Fig.  (3)  justifies  the  use  of 
the  scaling  law  based  on  the  Z  gradient  for  modeling  the  SGS  scalar  variance  (see  Section 
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VIC 2)  since  the  behavior  of  xt  and  \f  is  similar. 

Considering  the  r.m.s.  activity,  in  concert  with  the  findings  from  analysis  of  Eq.  (4), 
the  effect  of  the  Da />  variation  at  small  scales  cannot  be  underestimated.  For  HN600  and 
OHeGOO  the  SGS  turbulent  fluxes  are  larger  in  magnitude  than  the  molecular  diffusion  flux, 
meaning  that  the  scalar  fluctuations  contribute  to  the  mixing  at  the  SGS  scales.  Among  the 
four  SGS  diffusivity  terms,  the  second  one  assumes  large  values  when  compared  to  the  other 
diffusion-like  terms  or  advection  term,  indicating  that  it  must  be  retained  and  modeled  in 
the  SGS  scalar  variance  equation;  the  third  and  fourth  SGS  diffusivity  terms  have  smaller 
magnitudes  and  their  nature  is  dissipative;  finally,  the  first  one  compares  in  magnitude  to 
the  resolved  diffusion  in  HN600  and  OHeGOO  but  is  much  larger  than  the  resolved  diffusion 
for  OH750.  The  conclusion  is  that  particular  attention  should  be  devoted  to  SGS  Diffusivity 
1  and  SGS  Diffusivity  2  terms  because  they  have  a  dominant  contribution  to  the  mixing  of 
oxygen  with  hydrogen  at  the  SGS  scales.  The  Soret  contribution,  labeled  SGS  Soret  1,  has 
a  diffusion- like  behavior  and  is  non  negligible  only  in  the  HN600  mixing  layer.  Also,  the 
SGS  contributions  stemming  from  filtering  of  the  pressure-gradient  terms  are  negligible  in 
all  cases. 

B.  Assessment  of  the  presumed-PDF  approach  for  computing  filtered  non-linear 
scalar- dependent  quantities 

An  a  priori  evaluation  of  the  statistical  SGS  PDF  for  the  HN600  layer  is  illustrated  in 
Fig.  (4)  for  the  same  two  (xi,xs)  homogeneous  planes  of  coordinates  X2/Su,o  =  0.44  and 
^2/^,0  =  5.11  as  in  Fig.  (3).  The  DNS-extracted  filtered  PDF  <  fsgSc  >c  is  evaluated  by 
mass-weighted  averaging  the  Z  PDF  over  coarsened-grid  planes.  The  presumed  averaged 
SGS  PDFs  are  computed  employing  mapped  LES  moments  (£,  og)  and  using  Eq.  (15).  As 
expected,  the  figure  displays  the  best  prediction  when  a  /3  PDF  is  used.  The  difference 
between  one-moment,  and  two-moment  distributions  is  evident  both  in  the  central  part  of 
the  mixing  layer  (a;2/cG,o  =  0.44)  and  away  from  it  =  5.11).  The  difference  between 

the  two-moment,  distributions,  Gaussian  and  3.  is  also  clear,  particularly  in  zones  of  poor 
mixing  where  a  Gaussian  distribution  is  not  appropriate. 

Figure  (5)  portrays  the  corresponding  OH750  results  over  the  same  x-i  planes  as  for  the 
HN600  layer.  In  the  OH750  case,  no  large  difference  between  PDF  models  are  visible  which 
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is  conjectured  to  result  from  the  high  hydrogen  diffusivity  which  promotes  good  mixing. 
Finally,  OHeGOO  layer  results  are  shown  in  Figure  (6).  For  OHe  mixing,  the  two- moment 
PDFs  increase  the  model  accuracy  compared  to  the  single- moment  one,  and  a  B  PDF  gives 
slightly  better  results  than  the  Gaussian  PDF,  as  expected. 

Results  from  the  optimal  estimators  are  displayed  in  Fig.  (7).  Because  this  methodology 
is  based  on  conditioning  over  a  couple  of  moments,  the  Dirac  presumed  PDF  is  not  consid¬ 
ered.  In  order  to  increase  the  sample  size,  the  conditional  expectations  are  computed  over  a 
slab  of  the  mixing  layer  of  coordinates  x2/SUJt 0  =  ±5.3.  As  expected,  in  regions  of  well- mixed 
species,  the  (3  PDF  and  the  Gaussian  PDF  give  similar  predictions  while  the  /3  PDF  is  more 
appropriate  for  unmixed  situations. 

Thus,  we  have  shown  that  two  methodologies  are  in  agreement  regarding  the  appropri¬ 
ateness  of  the  presumed  (3  PDF  shape.  However,  it  is  known  that  for  atmospheric  conditions 
[17]  good  agreement  with  a  template  for  averaged  distributions  (as  given  by  both  method¬ 
ologies)  does  not  necessarily  imply  similar  pointwise  agreement,  so  it  is  important  to  test 
the  capabilities  of  the  three  distributions  of  Eqs.  (16),  (18)  and  (19)  to  provide  local  agree¬ 
ment  (as  needed  in  reactive  flows)  in  computations  where  they  would  be  used  to  reproduce 
filtered  non-linear  terms  (e.g.  reaction  rates).  For  this  purpose,  we  select  a  simple  non-linear 
Z-function 

F(0  =  exp  {-2  [erf-1  (2f  -  l)f}  ,  (50) 

representing  the  %  functional  form  in  a  one-dimensional  unsteady  laminar  subcritical  mixing 
layer  [47].  Then,  the  filtered  value  of  F(£)  may  be  estimated  by  integration  over  the  scalar 
space  using  Eq.  (13)  as 

bS=  f (6i) 

Jo 

- — mod 

Figure  (8)  illustrates  scatter  plots  of  the  modeled  quantity  F(Z )  (superscript  mod  labels  a 
modeled  quantity)  versus  the  exact  quantity  F(Z)  extracted  from  the  filtered  DNS  database 
HN600,  over  two  planes,  in  the  central  part  and  at  the  periphery  of  the  layer.  The  modeled 
functions  are  computed  by  replacing  fsgSc  (G  xj  in  Eq.  (51)  with  the  presumed  SGS  PDF 
(Dirac,  Gaussian  or  f3 )  constructed  from  the  exact  moments  extracted  from  the  filtered  DNS 
database  at  t*r.  The  local  results  are  in  accord  with  the  assessment  of  the  statistical  SGS 
PDFs  presented  in  Figs.  (4)  and  (7).  A  drastic  improvement  in  predictions  is  obtained  by 
using  the  /3-PDF  in  the  mixing  of  heptane  and  nitrogen.  An  improvement  is  also  obtained 
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for  the  oxygen/helium  system  (not  shown)  while  in  the  oxygen/hydrogen  mixing  layer  the 
differences  are  minor  (not  shown). 

As  mentioned  in  Section  III,  a  clipped  Gaussian  was  also  considered.  For  the  filtered 
PDF,  results  using  the  clipped  Gaussian  exhibited  small  differences  from  those  obtained 
with  the  truncated  Gaussian  for  the  HN600  mixing  layer  (not  shown),  and  imperceptible 
differences  for  the  OH750  and  OHeGOO  mixing  layers  (not  shown).  The  clipped  Gaussian 
also  provided  slightly  more  accurate  results  when  used  to  model  the  non-linear  Z-function 
(not  shown). 

Whereas  the  Z-function  of  Eq.  (50)  is  here  used  only  to  test  the  ability  of  the  models  to 
reproduce  nonlinearities,  the  function  is  not  necessarily  expected  to  represent  the  real  y  for 
three-dimensional  supercritical  flows.  This  topic  is  addressed  in  Appendix  A. 

C.  Evaluation  of  direct  models  for  the  SGS  scalar  variance 

Having  shown  in  Section  VIA  that  under  supercritical  p  conditions  modeling  of  new 
terms  in  the  az  equation  is  necessary  (but  uncertain),  and  shown  in  Section  VI B  that  well- 
known  assumed  f3  PDFs  may  be  used  in  modeling  nonlinear  scalar-dependent-  functions,  the 
next  step  is  to  assess  models  for  computing  az  from  the  filtered  DNS  solution  to  enable 
the  construction  of  the  presumed  SGS  PDFs.  Results  from  two  such  models  are  presented 
below. 

1.  Evaluation  of  the  approximate  deconvolution  model 

All  presented  results  were  computed  on  the  DNS  (rather  than  LES)  grid.  In  all  figures, 
the  information  is  shown  for  several  A/  Ax  dns  values  and  for  each  A/  Ax  dns  value  for  five 
orders  of  reconstruction. 

Since  Eq.  (24)  is  at  the  core  of  the  compressible  ADM  procedure,  we  first  inquired  about 
the  convergence  of  Z**  to  Z  according  to  the  approximation  of  Eq.  (25).  To  explore  this 
issue,  we  first  computed  the  local  second-order  moment  of  Z*  and  compared  it  to  the  exact 
moment  extracted  from  the  filtered  DNS  database  of  the  HN600  mixing  layer  at  t*r ;  the 
results  are  illustrated  in  Fig.  (9).  Then,  we  calculated  the  local  second-order  moment  of 
Z**  and  compared  it  to  the  exact  moment  extracted  from  the  filtered  DNS  database  of  the 
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HN600  mixing  layer  at  t,*r\  the  results  are  illustrated  in  Fig.  (10).  Whereas  the  moment 
of  Z*  converges  to  that  of  Z  as  N  increases,  the  observations  are  that  the  moment  of  the 
deconvoluted  approximate  field  Z**  does  not  converge  to  that  of  Z.  and  that  in  some  zones 
this  second-order  moment  is  overestimated.  This  overestimate  occurs  for  a  filter  width  as 
small  as  6A xdns  and  the  situation  deteriorates  with  increasing  values  of  A/A xdns-  The 
r.m.s.  activity  of  the  difference  between  Z**  and  Z*  is  presented  in  Fig.  (11).  Clearly,  this 
difference  increases  with  the  filter  width  and  it  also  globally  increases  with  N.  In  order 
to  inquire  about  the  impact  of  the  Z**  ~  Z*  approximation  on  the  SGS-scalar- variance 
predictions  for  compressible  flows,  the  ADM  procedure  is  tested  using  the  Favre  SGS  scalar 
variance. 

The  results  of  Fig.  (12)  portray  the  ADM  procedure  applied  to  the  primitive  variable  Z 
as  in  Eq.  (21),  and  the  SGS  scalar  variance  is  computed  using  the  exact  local  density  and 
the  approximate  field  Z*  as  az  =  {pZ*Z* /~p)  —  ( pZ*  /~p )  (pZ* /jo).  The  results  of  Fig.  (13) 
differ  from  those  of  Fig.  (12)  in  that  the  ADM  procedure  is  applied  to  both  conservative 
variables  p  and  pZ,  as  one  would  do  in  a  practical  case,  and  the  reconstructed  density 
p*  with  the  approximate  field  Z**  are  used  to  compute  the  SGS  scalar  variance  according 
to  Eq.  (26).  Comparing  Figs.  (12)  and  (13)  at  same  A/Axdns  value,  it  is  clear  that 
using  the  ADM  on  the  conservative  variables  improves  the  model  predictions.  For  example, 
for  A/Axdns  =  8,  a  third-order  approximation  is  sufficient  for  reproducing  (p(TZ)  at  the 
center  of  the  mixing  layer,  while  by  using  the  primitive  quantities  only  90%  of  its  value  is 
recovered.  The  difference  between  the  ADM  based  on  primitive  variables  and  that-  based  on 
conservative  variables  becomes  enlarged  as  the  filter  width  increases.  For  small  filter  widths 
(A/A xdns  =  2,4,6),  a  second-order  reconstruction  (3  terms)  gives  very  good  agreement- 
compared  t-o  the  exact  value.  Similarly  t-o  the  observation  of  Pant-ano  and  Sarkar  [16]  for 
incompressible  flows,  the  ADM  accuracy  depends  at-  least  on  the  level  of  turbulence,  i.e. 
the  Reynolds  number,  and  on  the  filter  width;  here,  there  is  the  additional  complication  of 
variable  density  which  introduces  Eq.  (24).  Figures  (14)  and  (15)  show  corresponding  results 
for  OH750  and  OHe600,  respectively,  when  the  ADM  is  performed  using  the  conservative 
variables. 

Despite  the  better  performance  of  the  ADM  conservative- variable  based  model  compared 
t-o  the  primitive-variable  based  one  for  the  third-order  approximation,  the  convergence  issue 
discussed  above  is  still  an  item  of  concern  when  using  the  approximation  Eq.  (25).  One 


24 


conclusion  from  the  presented  assessment  is  that  users  of  the  ADM  should  be  cautious 
when  employing  this  methodology  for  conservative  variables  in  the  context  of  compressible 
flows,  and  results  should  be  carefully  verified.  Here,  because  the  overestimation  of  the  SGS 
second-order  moment  of  Z  is  combined  with  an  underestimation  of  p.  the  global  result  is  a 
satisfactory  prediction  of  the  SGS  (Favre)  scalar  variance.  However,  ADM  should  only  be 
considered  here  as  an  approximation  rather  than  an  asymptotically  convergent  expansion. 

2.  Evaluation  of  the  dynamic  gradient-based  model 

Figure  (16)  depicts  planar  averages  of  the  modeled  SGS  scalar  variance  conditioned  on  the 
exact  one  extracted  from  the  filtered  HN600  DNS  database.  The  plots  represent  averages  at 
t*tr  over  a  plane  close  to  the  center  of  the  mixing  layer  o  =  0.44)  as  a  function  of  J)az, 

for  several  filter  widths;  the  vertical  arrow  represents  (paz)  and  provides  an  indication  of 
the  model  fidelity  at  that  particular  value.  For  A/ Axdns  =  2,  filtering  is  clearly  performed 
in  the  dissipation  range,  whereas,  as  an  example,  for  A/ Axdns  =  14,16  test-filtering  is 
performed  close  to  the  production  range,  and  thus  neither  of  these  values  are  in  concert 
with  SGS  modeling  assumptions,  but  they  are  here  presented  for  illustrative  purposes.  At 
A/ Axdns  =  2,  the  model  based  on  the  Leonard  term  expansion  agrees  with  the  classical 
model  used  in  conjunction  with  the  correct  filter  width,  but  neither  one  of  the  models  repro¬ 
duces  the  exact  value,  which  is  better  rendered  by  the  classical  model  utilized  in  conjunction 
with  an  incorrect  filter  width;  this  result  should  serve  as  a  warning  that  if  SGS  modeling  is 
tested  in  the  incorrect  wavenumber  range,  results  from  this  test  are  not  necessarily  reliable. 
Over  the  A/  Axdns  =  4  to  8  range,  the  new  model  and  the  classical  model  using  the  correct 
filter  width  agree  and  additionally  reproduce  the  DNS-extracted  value,  whereas  the  classical 
model  utilized  with  the  incorrect  filter  width  overpredicts  it;  however,  the  (p&z)  value  is 
equally  well  predicted  by  all  models.  For  A/ Axdns  =  10, 12,  the  fidelity  of  the  new  model 
to  predict  the  exact  is  maintained,  but  that  of  the  classical  model  with  the  correct 

filter  deteriorates  by  under  predicting  the  template,  and  for  A/  Axdns  =  12  severe  underpre¬ 
dictions  are  obtained  with  the  classical  model  in  conjunction  with  the  incorrect  filter  width. 
Most  important,  the  classical  model  used  with  either  filter  widths  produces  incorrect  values 
even  for  (paz),  whereas  the  new  model  maintains  high  fidelity  for  this  quantity.  The  robust¬ 
ness  of  the  new  model  is  highlighted  by  the  A/ Axdns  =  14, 16  results  where  its  predictions 
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are  still  excellent  whereas  those  with  the  classical  model  used  with  the  correct  filter  width 
display  severe  deteriorations  and  those  of  the  classical  model  using  the  incorrect  filter  width 
are  totally  compromised  by  producing  negative  values  of  the  SGS  scalar  variance.  Similar 
to  the  A/ Ax dns  =  10,12  situation,  {paz)  is  correctly  only  predicted  by  the  new  model, 
and  negative  values  of  the  SGS  scalar  variance  are  exhibited  by  the  classical  model  in  con¬ 
junction  with  the  incorrect  filter  width.  To  show  that  these  comparisons  are  not  .ly-plane 
dependent,  scatter  plots  at  t*tr  of  modeled  SGS  scalar  variances  against  actual  values  for 
several  filter  widths,  over  a  plane  close  to  the  periphery  of  the  mixing  layer  (£2/^,0  =  5.11) 
are  illustrated  in  Fig.  (17).  Only  for  A/A xdns  =  8  do  the  scatter  plots  overlap  for  all  three 
models,  and  predicted  negative  variances  by  the  classical  model  using  the  incorrect  filter 
width  appear  for  a  filter  ratio  of  A/  Ax  dns  =  12  which  is  smaller  than  the  14  and  16  where 
we  found  negative  values  in  the  =  0.44  plane. 

The  Figs.  (16)  and  (17)  results  were  for  two  selected  x2  planes.  To  further  evaluate  the 
potential  of  the  various  models  we  illustrate  in  Fig.  (18)  comparisons  of  the  modeled  (p(TZ) 
with  the  exact  one  for  the  same  A/ Ax  dns  values  and  for  the  entire  x2  significant  range. 
The  advantage  of  the  new  model  over  the  classical  model  using  the  correct  filter  width  is 
evident  for  as  small  value  as  A/  Ax  dns  =  8,  and  comparing  with  the  classical  model  using 
the  incorrect  filter  width  for  as  small  value  as  A/ Ax  dns  =  4.  The  high  fidelity  of  the  new 
model  persists  at  large  A/ Ax  dns  whereas  it  substantially  deteriorates  for  the  other  two 
models  with  increasing  A/ Ax  dns  values.  The  model  coefficients  computed  with  the  three 
models  are  depicted  in  Fig.  (19).  The  indications  are  that  over  all  ^2-planes  of  the  mixing 
layer,  the  use  of  the  Leonard  term  expansion  for  the  dynamic  model  yields  model  coefficient 
values  which  span  a  smaller  range  at  fixed  A/  Ax  dns  value  than  those  of  the  classical  model, 
thereby  showing  greater  potential  in  maintaining  stability  of  a  LES  computation. 

Corresponding  results  for  the  OH750  database  are  displayed  in  Fig.  (20),  (21),  (22),  and 
for  the  OHe600  database  are  depicted  in  Figs.  (23),  (24),  (25).  For  the  OH750  conditional 
averages  and  mean  profiles,  the  advantage  of  the  new  model  predictions  are  less  drastic 
when  compared  to  the  classical  model  than  for  the  HN600  database.  The  new  model  and 
the  classical  one  with  the  correct  filter  width  are  in  close  agreement,  while  their  results 
differ  from  those  using  the  classical  procedure  with  the  incorrect  filter  width  A  =  A.  For  the 
OHe600  database  which  has  the  highest  Retr  value  among  the  three  examined,  the  results 
at  small  to  moderate  A/ Ax  dns  values  are  the  same  as  for  HN600  for  the  o  =  0.44 
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conditional  averages,  however,  for  large  A/  Axdns  values  close  to  the  production  range  the 
predictions  from  all  three  models  at  this  particular  ^-coordinate  converge  and  underestimate 
the  exact  value.  The  advantage  of  the  new  model  is  though  evident  for  predicting  the  mean 
SGS  variance  value  (Fig.  (24))  over  all  x2-planes  of  the  mixing  layer  even  at  large  A/A x^ns 
values. 

D.  Modeled  scalar  variance  in  conjunction  with  the  presumed-PDF  approach 

As  a  recall,  the  results  of  Fig.  (8)  were  obtained  with  the  exact  moments  of  the  SGS  PDF, 
as  extracted  from  the  filtered  DNS  database.  It  is  thus  important  to  explore  the  potential 
of  the  two  SGS-scalar- variance  models  examined  in  Section  VI C,  namely  the  ADM  using  a 
third-order  approximation  and  the  new  dynamic  gradient-like  model.  To  this  end,  we  use 
each  of  these  models  to  construct  the  3-  PDF  and  assess  their  performance  in  reproducing 
the  same  filtered  non-linear  A-function  of  Eq.  (50),  as  in  Section  VI B.  Predictions  are 
illustrated  in  Fig.  (26)  for  HN600  at  A/A xdns  =  8  and  £2/^,0  =  0.44,  as  an  example. 
Each  of  the  models  is  evaluated  through  scatter  plots.  Independent  of  the  SGS-scalar- 
variance  model  used,  the  predictions  are  excellent,  showing  a  priori  the  potential  of  the  new 
(Jz  models  combined  with  a  PDF  approach  to  reproduce  filtered  non-linear  functions  of  the 
scalar,  as  would  be  the  dissipation  rate,  reaction  rates,  etc.  Examining  the  scatter  plots,  the 
ADM  results  show  a  smaller  dispersion  than  the  gradient-based  model.  On  the  other  hand, 
the  former  is  more  expensive  than  the  latter  and  the  ADM  exhibits  convergence  problems 
of  the  deconvolution  series  for  compressible  flow,  casting  uncertainty  on  success  at  high 
Reynolds  numbers  where  convergence  will  be  further  influenced  by  turbulence-dependent 
aspects. 

VII.  SUMMARY  AND  CONCLUSIONS 

The  goal  of  this  study  was  to  investigate  the  modeling  of  the  SGS  scalar  variance  un¬ 
der  supercritical-pressure  conditions  where  the  real-gas  equation  of  state,  the  full  (3-term) 
expression  for  the  species  mass  diffusion  flux  and  transport  properties  varying  with  the  ther¬ 
modynamic  variables,  all  preclude  assuming  that  the  same  models  as  those  at  atmospheric- 
pressure  conditions  are  valid.  To  this  end,  we  followed  the  classical  approach  whereby 
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the  SGS  scalar  variance  equation  is  derived  and  terms  non-computable  from  the  LES  so¬ 
lution  are  modeled,  with  the  intent  of  providing  closure  and  solving  the  equation.  Thus, 
we  first  developed  the  equation  describing  the  evolution  of  the  SGS  scalar  variance  un¬ 
der  supercritical-pressure  conditions  and  highlighted  its  complexity,  particularly  recognizing 
terms  not  present  under  atmospheric-pressure  conditions.  We  also  presented  a  second  form 
of  the  equation  which  was  more  adept  at  highlighting  the  nature  of  the  new  SGS  terms  and 
their  contributions.  The  activity  of  terms  in  the  first  form  of  the  equation  was  examined 
using  a  filtered  DNS  database  represented  by  transitional  states  describing  the  mixing  of 
binary  species,  for  three  systems  of  species,  under  supercritical  pressure  conditions.  The 
findings  were  that  the  activity  of  some  of  these  new  terms  is  of  same  magnitude  as  that 
of  classical  terms,  meaning  that  they  cannot  be  neglected.  Most  important  among  these 
new  terms  were  those  expressing  subgrid  activity  due  to  spatially  variable  diffusion  coeffi¬ 
cients.  The  second  form  of  the  equation  confirmed  the  importance  of  the  SGS  diffusivity 
and  identified  for  this  equation  a  new  dissipation  contribution  arising  from  the  Soret  term. 
Recognizing  that  no  SGS  models  are  available  to  model  these  terms,  and  thus  to  close  the 
SGS  scalar  variance  equation,  the  attention  was  instead  refocused  on  a  second  method,  that 
is,  the  direct  modeling  of  the  SGS  scalar  variance. 

This  second  route  first  involved  examining  the  SGS  PDF  of  the  scalar  by  assuming  the 
form  of  the  PDF  and  using  the  same  filtered  DNS  database  to  extract  the  exact  moments  of 
the  PDF.  Three  PDF  forms  were  investigated  -  the  Dirac,  Gaussian  and  (3  PDF  -  and  the 
results  showed  that  they  ranked  in  increasing  success  in  the  order  cited.  This  encouraging 
ability  of  the  3  PDF  motivated  the  development  of  two  direct  models  for  the  SGS  scalar  vari¬ 
ance.  The  first  SGS-scalar-variance  model  was  based  on  the  ADM  procedure  reformulated 
for  application  to  compressible  flows.  The  second  SGS-scalar-variance  model  was  based  on 
a  gradient-like  dynamic  model  using  the  Leonard  term  expansion.  Success  with  these  two 
models  motivated  a  reassessment  of  the  ability  to  model  a  filtered  non-linear  function  of 
the  scalar  by  using  the  3  PDF  in  conjunction  with  either  one  of  these  models  for  the  SGS 
scalar  variance,  and  with  the  mean  computed  from  the  filtered  DNS.  The  findings  were  that 
either  one  of  the  direct  SGS  scalar  variance  models  provided  a  high-fidelity  duplication  of 
the  DNS-extracted  SGS  PDF,  which  is  manifested  by  the  excellent  reproduction  of  a  filtered 
non-linear  function  of  the  scalar.  Although  the  ADM  was  generally  more  accurate  than  the 
gradient-based  model,  it  was  shown  that  the  ADM  procedure  is  not  necessarily  convergent 


28 


for  compressible  flows;  thus,  the  results  could  be  interpreted  as  an  approximation  rather 
than  a  model  the  accuracy  of  which  asymptotically  increases  with  series  higher  truncation 
order. 

Therefore,  although  supercritical-pressure  conditions  entail  new  challenges  in  the  model¬ 
ing  of  the  SGS  scalar  variance,  these  challenges  were  met  for  describing  the  mixing  of  several 
binary- species  systems.  Further  a  posteriori  studies  should  reveal  the  true  potential  of  these 
models,  and  applications  to  reacting  flows  would  represent  the  ultimate  test. 


APPENDIX  A:  ASSESSMENT  OF  THE  ONE-DIMENSIONAL  LAMINAR 
SCALAR  DISSIPATION  FOR  TRANSITIONAL  SUPERCRITICAL-PRESSURE 
MIXING  LAYERS 

The  dissipation  rate  in  the  context  of  the  flamelet  model  is  often  modeled  in  the  framework 
of  a  one-dimensional  counterflow  [48].  Under  this  assumption,  an  analytical  solution  for  the 
dissipation  rate  is  available  under  subcritical  conditions  [1]  which  is  the  functional  form 
of  Eq.  (50).  In  order  to  evaluate  whether  this  form  still  holds  in  the  present  case,  the 
averages  of  both  \t  and  xf  of  Eqs.  (48)  and  (49)  are  computed  from  the  DNS  conditioned 
on  F(Z)  (Eq.  50).  The  results  are  illustrated  in  Fig.  (27).  For  both  \T  and  xf,  there  is 
a  linear  dependency  on  F(Z)  in  zones  corresponding  to  large  values  of  Z.  while  departures 
are  observed  for  smaller  values  of  Z.  This  indicates  that  F(Z)  is  unreliable  since  it  cannot- 
handle  regions  of  fluid  mixing  where  the  gradients  would  be  largest. 
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FIG.  1:  Planar  r.m.s.  activity  of  terms  in  Eq.  (4)  as  extracted  from  the  filtered  DNS  databases  at 
tfr.  Top:  HN600;  center:  OH750;  bottom:  OHe600.  A/Axdns  =  8.  Units  are  kg/(m3  s). 
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FIG.  2:  Planar  averages  (left)  and  planar  r.rn.s.  activity  (right)  of  terms  in  Eq.  (9).  Extracted  from 
the  filtered  DNS  databases  at  t*tr.  Top:  HN600;  center:  OH750;  bottom:  OHe600.  A/A xdns  =  8. 
Units  are  kg/(m3  s). 
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<  Xt\xf  > 


FIG.  3:  Conditional  average  of  xt  on  xf  extracted  from  the  HN600  DNS  database  over  planes  of 
coordinates  X2/<L,o  =  0.44  (left)  and  x' 2/^,0  =  5.11  (right)  at  tfr.  A/Axdns  =  8. 
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FIG.  4:  Filtered  PDFs  <  fSgsc  >o  from  the  HN600  DNS  at  over  planes  of  coordinates  X2/Su>,0  = 
0.44  (left)  and  £2/^,0  =  5.11  (right).  Exact  SGS:  solid  lines.  Assumed  PDFs  are  the  Dirac 
(top),  the  Gaussian  (center)  and  the  Beta  (bottom)  functions.  A/Axdns  =  4,  dashed  lines; 
A/Axdns  =  8,  dotted  lines;  A/Axdns  =  12,  dot-dashed  lines. 
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FIG.  5:  Filtered  PDFs  <  fsgSc  >c,  from  the  OH750  database  at  over  planes  of  coordinates 
*2/^,0  =  0.44  (left)  and  X2/8ufi  =  5.11  (right).  Exact  SGS:  solid  lines.  Assumed  PDFs  are  the 
Dirac  (top),  the  Gaussian  (center)  and  the  Beta  (bottom)  functions.  A/Axdns  =  4,  dashed  lines; 
A/Axdns  =  8,  dotted  lines;  A/Axdns  =  12,  dot-dashed  lines. 
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FIG.  6:  Filtered  PDFs  <  fsgSc  >c,  from  the  OHe600  database  at  t^r  over  planes  of  coordinates 
*2/^,0  =  0-44  (left)  and  X2/8ufi  =  5.11  (right).  Exact  SGS:  solid  lines.  Assumed  PDFs  are  the 
Dirac  (top),  the  Gaussian  (center)  and  the  Beta  (bottom)  functions.  A/Axdns  =  4,  dashed  lines; 
A/Axdns  =  8,  dotted  lines;  A/Axdns  =  12,  dot-dashed  lines. 
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FIG.  7:  Conditional  expectation  of  fsgi 5c(£;£,  cr^)  over  (£,  o^)  and  over  the  slab  at  £2/^,0  =  ±5.3. 
Top:  HN600,  £  =  0.805  ±  0.005,  =  0.0505  ±  0.0005  (left)  and  £  =  0.505  ±  0.005,  =  0.0105  ± 

0.0005  (right).  Center:  OH750,  £  =  0.905  ±  0.005,  ct5  =  0.00505  ±  0.00005  (left)  and  £  =  0.505  ± 
0.005,  cre  =  0.00505  ±  0.00005.  (right).  Bottom:  OHe600,  £  =  0.905  ±  0.005,  =  0.0105  ±  0.0005 
(left)  and  £  =  0.505  ±  0.005,  =  0.0105  ±  0.0Q£jj  (right).  A/ A xdns  =  8. 


F(Z) 
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FIG.  8:  Scatter  plot  of  the  modeled  function  F(Z)  versus  the  exact  quantity  F(Z)  computed 
over  planes  X2/5u,o  =  0.44  (left)  and  £2/^,0  =  5.11  (right).  The  exact  quantity  is  the  filtered 
HN600  DNS  at  t*r  and  the  models  are  the  Dirac  PDF  (top),  Gaussian  PDF  (center)  and  f3  PDF 
(bottom).  A/ Ax dns  =  8. 
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%2/  <L,0 


X2/ Su,0 


FIG.  9:  Predictions  of  the  SGS  scalar  variance  using  the  deconvoluted  field  Z*.  Several  orders  of 
approximation  are  shown  for  different  filter  widths  at  t^r  for  the  HN600  mixing  layer.  Exact  SGS 
scalar  variance:  solid  line.  Dotted  lines  are  the  first  to  fifth  order  approximations,  the  third  order 
being  distinguished  by  a  dash-dotted  line.  Variances  are  non-dimensionalized  by  the  exact  value 
at  the  center  of  the  mixing  layer. 
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FIG.  10:  Predictions  of  the  SGS  scalar  variance  using  the  deconvoluted  field  Z**.  Several  orders  of 
approximation  are  shown  for  different  filter  widths  at  t^r  for  the  HN600  mixing  layer.  Exact  SGS 
scalar  variance:  solid  line.  Dotted  lines  are  the  first  to  fifth  order  approximations,  the  third  order 
being  distinguished  by  a  dash-dotted  line.  Variances  are  non-dimensionalized  by  the  exact  value 
at  the  center  of  the  mixing  layer. 
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FIG.  11:  Planar  r.m.s.  activity  of  the  field  ( Z **  —  Z*)  for  several  orders  of  ADM  approximation 
and  different  filter  widths  at  t*r  for  the  HN600  mixing  layer.  Thin  solid  line:  first  order  of  approx¬ 
imation;  thick  solid  line:  fifth  order  of  approximation.  Dotted  lines  are  the  second  to  fourth  order 
approximations. 
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FIG.  12:  Predictions  of  the  Favre  SGS  scalar  variance  using  ADM  applied  to  the  primitive  quantity. 
Several  orders  of  approximation  are  shown  for  different  filter  widths  at  t*r  for  the  HN600  mixing 
layer.  Exact  SGS  scalar  variance:  solid  line.  Dotted  lines  are  the  first  to  fifth  order  approximations, 
the  third  order  being  distinguished  by  a  dash-dotted  line.  Variances  are  non-dimensionalized  by 
the  exact  value  at  the  center  of  the  mixing  layeiL 


FIG.  13:  Predictions  of  the  Favre  SGS  scalar  variance  using  ADM  applied  to  conservative  quan¬ 
tities.  Several  orders  of  approximation  are  shown  for  different  filter  widths  at  t£r  for  the  HN600 
mixing  layer.  Exact  SGS  scalar  variance:  solid  line.  Dotted  lines  are  the  first  to  fifth  order 
approximations,  the  third  order  being  distinguished  by  a  dash-dotted  line.  Variances  are  non- 
dimensionalized  by  the  exact  value  at  the  center,  of  the  mixing  layer. 
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FIG.  14:  Predictions  of  the  Favre  SGS  scalar  variance  using  ADM  applied  to  conservative  quan¬ 
tities.  Several  orders  of  approximation  are  shown  for  different  filter  widths  at  for  the  OH750 
mixing  layer.  Exact  SGS  scalar  variance:  solid  line.  Dotted  lines  are  the  first  to  fifth  order 
approximations,  the  third  order  being  distinguished  by  a  dash-dotted  line.  Variances  are  non- 
dimensionalized  by  the  exact  value  at  the  centejLof  the  mixing  layer. 


FIG.  15:  Predictions  of  the  Favre  SGS  scalar  variance  using  ADM  applied  to  conservative  quan¬ 
tities.  Several  orders  of  approximation  are  shown  for  different  filter  widths  at  t\r  for  the  OHe600 
mixing  layer.  Exact  SGS  scalar  variance:  solid  line.  Dotted  lines  are  the  first  to  fifth  order 
approximations,  the  third  order  being  distinguished  by  a  dash-dotted  line.  Variances  are  non- 
dimensionalized  by  the  exact  value  at  the  center,  of  the  mixing  layer. 
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FIG.  16:  Planar  averages  of  modeled  SGS  scalar  variances  conditioned  on  the  exact  ones,  evaluated 

using  the  filtered  HN600  DNS  at  tjfr,  over  the  plane  £2/^,0  =  0.44.  Coefficients  are  computed  using 

— 2  -  —2 

the  new  model  (Eq.(47))  (squares),  and  the  classical  model  (Eq.(46))  with  either  A  =  A2  +  A 
(empty  circles)  or  A  =  A  (filled  circles).  The  arrow  indicates  the  mean  quantity  ( ~paz )  ■ 
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FIG.  17:  Scatter  plot  of  modeled  against  exact  SGS  scalar  variance  evaluated  using  the  filtered 

HN600  DNS  at  t*r,  over  the  plane  X2/SUjo  =  5.11.  Coefficients  are  computed  using  the  new  model 

— 2  -  —2 

(Eq.(47))  (blue  symbols),  and  the  classical  model  (Eq.(46))  with  either  A  =  A2  +  A  (black 
symbols)  or  A  =  A  (red  symbols).  R  =  A/ A xdns- 
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FIG.  18:  Profiles  of  modeled  and  exact  SGS  scalar  variances,  evaluated  using  the  filtered  HN600 

DNS  at  t%r.  Exact  values:  solid  line.  New  model:  dash-dotted  line.  Classical  model  using 

— 2  — 2  — 

A  =  A2  +  A"  :  dotted  line.  Classical  model  using  A  =  A  :  dashed  line.  Variances  are  non- 
dimensionalized  by  the  exact  value  at  the  center  of  the  mixing  layer. 
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FIG.  19:  Model  coefficients  Cc[  evaluated  using  the  filtered  HN600  DNS  at  t£r.  New  model:  dash- 

— 2  ~  — 2  — 

dotted  line.  Classical  model  using  with  A  =  A2  +  A  :  dotted  line.  Classical  model  using  A  =  A  : 

dashed  line. 
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FIG.  20:  Planar  averages  of  modeled  SGS  scalar  variances  conditioned  on  the  exact  ones  evaluated 

using  the  filtered  OH750  DNS  at  t\r  over  the  plane  £2/^,0  =  0.44.  Coefficients  are  computed  using 

—2  -  —2 

the  new  model  (Eq.(47))  (squares)  and  the  classical  model  (Eq.(46))  with  either  A  =  A2  +  A 
(empty  circles)  or  A  =  A  (filled  circles).  The  arrow  indicates  the  mean  quantity  (pcrz)  • 
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FIG.  21:  Profiles  of  modeled  and  exact  SGS  scalar  variances  evaluated  using  the  filtered  OH750 

DNS  at  t%r.  Exact  values:  solid  line.  New  model:  dash-dotted  line.  Classical  model  using 

— 2  — 2  — 

A  =  A2  +  A"  :  dotted  line.  Classical  model  using  A  =  A  :  dashed  line.  Variances  are  non- 
dimensionalized  by  the  exact  value  at  the  center  of  the  mixing  layer. 
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FIG.  22:  Model  coefficients  Cd  evaluated  using  the  filtered  OH750  DNS  at  t*r.  New  model:  dash- 

—2-—2  n 

dotted  line.  Classical  model  using  A  =  A2  +  A  :  dotted  line.  Classical  model  using  A  =  A.: 

dashed  line. 
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FIG.  23:  Planar  averages  of  modeled  SGS  scalar  variances  conditional  to  exact  ones  evaluated 

using  the  OHe600  DNS  at  t*r  over  the  plane  £2/<5uj,o  =  0.44.  Coefficients  are  computed  using  the 

— 2  -  —2 

new  model  (Eq.(47))  (squares)  and  the  classical  model  (Eq.(46))  with  either  A  =  A2  + A"  (empty 
circles)  or  A  =  A  (filled  circles).  The  arrow  indicates  the  mean  quantity  (paz)  ■ 
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FIG.  24:  Profiles  of  modeled  and  exact  SGS  scalar  variances  evaluated  using  the  filtered  OHe600 
DNS  at  tfr.  Exact  values:  solid  line.  New  model:  dash-dotted  line.  Classical  model  using 

— 2  ^o—2  — 

A  =  A2  +  A"  :  dotted  line.  Classical  model  using  A  =  A  :  dashed  line.  Variances  are  non- 
dimensionalized  by  the  exact  value  at  the  center  of  the  mixing  layer. 
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FIG.  25:  Model  coefficients  Cd  evaluated  by  using  the  filtered  OHe600  DNS  at  t*r.  New  model: 

—2-—2  — 
dash-dotted  line.  Classical  model  using  A  =  A2  +  A  :  dotted  line.  Classical  model  using  A  =  A  : 

dashed  line. 
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FIG.  26:  Scatter  plot  of  the  modeled  function  F(Z)  versus  the  exact  quantity  F(Z )  computed 
over  the  plane  X2/<^,o  =  0.44.  The  evaluation  is  made  for  the  filtered  HN600  DNS  at  t*r,  using 
the  scalar  variance  az  modeled  employing  a  third-order  ADM  (left),  and  the  new  dynamic  model 
(right).  A/ Axdns  =  8. 
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Thermomechanical  Response  of  a  Reactive  Gas  to  Extremely  Rapid  Transient, 
Spatially  Distributed  Energy  Addition:  An  Asymptotic  Formulation 


D.  R.  Kassoy1,  K.  Wojchiechowski2 

Abstract 

An  asymptotic  formulation  developed  for  inert  gas  heating  ((Kassoy,  DR.,  2010,  J.  Eng.Math,  68,  249-262).  is 
extended  to  describe  the  thermomechanical  response  of  a  reactive  gas  to  a  localized,  exothermic  chemical  reaction. 
The  model  is  developed  for  a  subcritical,  perfect  gas  with  one-step,  high  activation  energy  exothermic  kinetics.  A 
finite  volume  of  gas  (the  hot  spot)  is  heated  initially  by  an  external  source  on  a  time-scale  short  compared  to  the 
acoustic  time  of  the  region,  with  the  objective  of  raising  the  local  temperature  sufficiently  to  ignite  a  robust  rapid 
reaction.  The  analysis  defines  physical  and  chemical  conditions  compatible  with  nearly  inertially -confined  heating. 
Numerical  solutions  to  the  describing  equations  show  that  a  spatially  distributed  reaction  wave  appears 
spontaneously  in  the  hottest  portion  of  the  hot  spot,  propagates  initially  through  the  relatively  slowly  moving  fluid  at 
a  supersonic  speed  (relative  to  the  hot  gas  speed  of  sound)  and  then  decelerates  significantly  and  steepens 
considerably  as  the  front  nears  the  much  colder  boundary  of  the  hot  spot.  The  configuration  evolves  toward  a 
discontinuous  front  separating  hot,  high-pressure,  burned  gas  on  one  side  from  cold,  low  pressure  reactant  on  the 
other  side.  During  the  relevant  heating  time  scale  the  entire  process  occurs  in  a  nearly  incompressible  medium, 
leading  to  an  ephemeral,  isolated,  burned  out,  hot,  high  pressure  spot  embedded  in  a  cold  unreacted,  lower  pressure 
gas.  The  large  pressure  gradient  at  the  front  induces  a  local  positive  radial  fluid  velocity.  Fluid  expelled  from  the 
hot  spot  boundary  acts  as  source  of  mechanical  disturbances  propagating  into  the  neighboring  cold  reactant.  The 
amplitude  of  those  disturbances  depends  upon  the  energy  addition  level  during  the  reactive  phase  of  the  hot  spot. 
The  formulation  also  identifies  conditions  compatible  with  a  fully  compressible  heat  addition  process,  characterized 
by  a  very  significant  internal  expansion  Mach  number,  likely  the  source  of  reactive  blast  wave  generation  in  the 
environmental  reactant. 


1  Professor  Emeritus,  Mechanical  Engineering  Department,  University  of  Colorado,  Boulder,  CO  80309- 
0427,  USA  david.kassoy@colorado.edu 

2  Graduate  Student,  Mathematics,  University  of  Colorado,  Denver,  CO 


1.  Introduction 


Oppenheim  and  Soloukhin  [1]  state  that  “Gasdynamics  of  explosions  is  best  defined  as  the 
science  dealing  with  the  interrelationship  between  energy  transfer  occurring  at  a  high  rate  in  a 
compressible  medium  and  the  concomitant  motion  set  up  in  the  medium.”  This  perceptive 
insight  challenges  the  modeler  to  relate  quantitatively  transient,  spatially  resolved  energy 
deposition  into  a  hot  spot  to  changes  in  the  thermodynamic  state,  and  the  induced  fluid  motion, 
as  well  as  to  predict  the  acoustic  and/or  gasdynamic  disturbances  generated  in  the  unheated 
gaseous  environment  beyond  the  hot  spot  boundary. 

An  initial  step  toward  developing  a  comprehensive  thermomechanical  theory  for  hot  spots 
has  been  developed  by  Kassoy  [2].  An  external  source  provides  spatially  distributed,  transient 
heat  addition  to  a  finite  volume  of  inert  gas  on  a  time  scale  short  compared  to  the  acoustic 
time  of  the  heated  region.  Asymptotic  methods  are  used  to  derive  reduced  equation  systems 
that  describe  hot  spot  evolution  resulting  from  a  wide  range  of  energy  addition.  The  results 
demonstrate  that  nearly-inertially  confined  heat  addition  occurs  within  the  hot  spot  if  the 
amplitude  of  energy  is  less  than  a  critical  value.  This  type  of  nearly  constant  volume  heating  is 
characterized  by  a  synchronized  increase  in  spatially  distributed  temperature  and  pressure,  very 
small  changes  in  density  and  a  tiny  internal  expansion  Mach  number.  Hot  gas  expelled  from  the 
volume  at  a  locally  small  Mach  number  (the  “piston”  effect  [3])  during  the  heating  process 
generates  only  acoustics  in  the  neighboring  cold  gas  (the  far-field).  Sufficiently  large  energy 
addition  leads  to  a  fully  compressible  heat  addition  process.  Much  larger  “piston”  Mach 
numbers  are  the  sources  of  far-field  shock  and  blast  wave  propagation.  Extension  of  the 
theoretical  framework  to  reactive  gases  and  other  time  scales  is  the  subject  of  the  present  work. 

Localized  hot  spots  are  common  in  a  wide  range  of  reactive  gasdynamic  processes,  most 
notably  detonation  initiation  [1],  One  class  of  hot  spot  arises  from  relatively  rapid,  localized 
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thermal  energy  deposition  from  an  external  source  into  a  volume  of  reactive  gas  (e.g.,  laser  or 
spark).  Another,  the  “reaction  center”  [1]  arises  spontaneously  during  the  evolution  to 
detonation.  The  modeling  and  experimental  literature  on  the  subject  is  extensive,  and  is  referred 
to  within  the  limited  number  of  references  cited  in  this  work. 

Sufficiently  fast  and  large  energy  addition  can  facilitate  a  direct  initiation.  Clarke,  Kassoy 
and  Riley  [4]  and  Clarke,  Kassoy  Maharzi,  Riley  and  Vasantha  [5]  model  detonation  initiation 
following  spatially-resolved  transient  energy  deposition  from  a  planar  boundary  via 
conduction  into  an  adjacent  reactive  gas  in  a  semi-infinte  domain.  The  Navier-Stokes  equations 
are  integrated  numerically  to  reveal  a  quantitative  time-history  of  spatially  resolved  reactive 
gasdynamic  processes  that  lead  to  planar  detonation  formation.  The  authors  recognize  that 
“...direct  initiation  of  detonation  requires  sufficient  power  input  to  first  of  all  generate  a  suitably 
strong  precursor  shockwave,  which  then  becomes  the  trigger  to  switch  on  vigorous  chemical 
activity  in  its  wake.  The  hall  mark  of  this  vigor  is  its  capacity  to  exploit  the  inertia  of  the  fluid 
by  raising  local  pressures  and  temperatures,  with  little  diminution  in  local  density;  the  pressure 
waves  so  formed  propagate  and  increase  precursor  shock  strength  which  therefore  lifts  the 
overall  density  levels,  as  well  as  those  of  pressure  and  temperature.  All  of  these  processes 
interlock  in  a  continuously  accelerated  sequence  that  progresses  towards  a  steady  state.... ZND 
detonation.” 

Mazaheri  [6]  and  Eckett,  Quirk  and  Shepherd  [7,8],  model  planar  and  spherical  detonation 
initiation,  respectively,  initiated  by  blast  waves  subsequent  to  instantaneous  deposition  of 
energy  at  a  plane  or  a  point.  Energy  deposition  criteria  are  used  to  distinguish  between 
sustained  and  failed  detonations.  Computational  modeling  results  demonstrate  that  blast  wave 
propagation  through  unreacted  gas  mixture  leads  to  the  formation  of  localized  regions  of  rapid 


3 


chemical  heat  release  (thermal  explosions)  characterized  by  relatively  high  temperature  and 
pressure,  similar  to  those  found  in  Refs.  9-11.  These  hot  spots  are  the  subsequent  sources  of 
compression  waves  that  may  run  up  to  and  strengthen  the  blast  wave  front  enough  to  generate 
and  sustain  a  classical  detonation  (shock  coupled  to  a  reaction  zone). 

Sileem,  Kassoy  and  Hayashi  [12],  (SKH),  Kassoy,  Keuhn,  Nabity  and  Clarke  [13,14], 
(KKNC),  Regele,  Kassoy  and  Vasilyev  [15]  model  the  reactive  gas  response  to  relatively  smaller 
spatially  distributed,  transient  energy  deposition  into  a  finite  target  volume.  They  describe 
the  sequence  of  reactive  gasdynamic  events  occurring  in  a  deflagration  to  detonation  transition 
(DDT),  including  the  initially  heated  volume  and  the  spontaneous  appearance  of  local  hot  spots 
far  from  the  original  energy  deposition.  Computational  results,  based  on  MacCormack 
numerical  methods  with  fixed  grids  [12-14]  and  the  Adaptive  Wavelet  Collocation  Method  [15], 
resolve  hot  spot  transients  facilitating  the  detonation  initiation  process. 

Gu  et  al.  [16]  use  computational  solutions  to  the  Euler  equations  with  multistep  kinetics 
relevant  to  FE-CO-Air  and  IE-Air  mixtures  to  identify  five  distinct  modes  of  reaction  front 
propagation  arising  from  a  preexisting  local  hot  spot.  They  find  that  evolution  of  the  detonative 
mode  depends  on  the  temperature  distribution  properties  of  the  hot  spot  and  “...the  ratio  of  the 
hotspot  acoustic  time  to  the  heat  release  rate  excitation  time. . .” 

Detonation  initiation  associated  with  reflected  shocks  and  shock  flame  interactions,  studied 
intensively  by  Oran  and  co-workers  beginning  in  the  1980’s,  and  reviewed  by  Oran  and  Gamezo 
[17]  is  also  affected  substantially  by  hot  spot  dynamics.  In  fact,  nearly  all  studies  of  detonation 
initiation  contain  qualitative  descriptions  of  the  role  played  by  hot  spots  in  the  development  of 
detonations.  Since  the  early  experimental  observation  by  Oppenheim  [18]  of  an  “explosion  in  an 
explosion”  (reaction  center)  it  has  been  argued  that  hot  spots  are  local  sources  of  compression 
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waves  that  strengthen  existing  lead  shocks  to  promote  the  existence  of  coupled  reaction  zones. 
That  argument  depends  on  the  simultaneous  local  increase  in  temperature  and  pressure  caused  by 
localized  chemical  heat  release  in  a  nearly  inertially  confined  fluid  volume  (constant  volume 
heating)  characterized  by  a  minimal  change  in  density. 

Zeldovich  and  co-workers  [19,20]  articulate  an  intuitive  theory  for  the  initiation  and 
“... (spontaneous)  propagation  of  intense  chemical  reaction...”  arising  from  a  pre-existing  hot 
spot  with  a  localized  linear  temperature  gradient.  The  authors  describe  a  model  for  constant 
volume  thennal  explosion  phenomena  (“...in  each  particle  of  the  substance,  thennal  explosion 
occurs  independently.”).  The  chemico-physical  conditions  necessary  for  this  type  of  reaction 
wave  initiation  and  propagation  are  not  articulated. 

Many  researchers  have  endeavored  to  rationalize  the  Zeldovich  model  concept.  Jackson, 
Kapila  and  Stewart  [9],  Short  [10,11],  Kapila,  Shwendeman,  Quirk  and  Hawa  [21],  have  used  a 
combination  of  asymptotic  and  computational  methodologies  to  model  the  evolution  of  spatially 
distributed  thermal  explosions  arising  from  a  pre-existing  localized  linear  temperature 
gradient  in  a  rigid  container.  Asymptotic  results  describe  spatially  distributed,  transient 
chemical  heat  release  during  a  relatively  lengthy  induction  period,  followed  by  much  shorter 
periods  of  extremely  rapid  energy  deposition  (the  explosion)  into  a  relatively  tiny  spatial  volume. 
Comprehensive  asymptotic  and  numerical  results  in  Ref.  21  provide  an  excellent  assessment  of 
the  impact  of  the  temperature  gradient  slope  on  the  combustion  wave  evolution.  Many  of  the 
observed  effects  can  be  attributed  to  differences  in  the  local  chemical  heat  release  time  and  the 
local  acoustic  time,  exploited  formally  in  Ref.  2.  These  confined  hot  spot  theories  [9-11,21] 
describe  the  thermomechanical  properties  of  the  chemically  heated  gas,  but  are  not  formulated  to 
describe  the  impact  of  focused  rapid  heat  release  on  gases  external  to  the  heated  region. 
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The  current  analysis  extends  the  asymptotic  methodology  in  Ref.  2  to  a  reactive  gas  with 
one  step,  high  activation  energy  Arrhenius  kinetics.  Initially,  an  external  source  provides 
thermal  power  to  a  finite  volume  of  gas  on  a  time  scale  short  compared  to  the  acoustic  time  of 
the  heat  region.  The  resulting  near-constant  volume  heat  addition  causes  a  relatively  modest 
temperature  increase,  insufficient  to  initiate  the  high  activation  energy  reaction.  Continued 
power  deposition  on  a  slightly  longer  time  scale,  but  shorter  than  the  acoustic  time  of  the 
volume,  leads  to  a  temperature  increase  large  enough  to  facilitate  a  strong  reaction.  This  strategy 
obviates  the  need  for  a  more  traditional  induction  period  thermal  explosion  analysis.  The 
reaction  evolution  and  subsequent  front  propagation  occurs  in  a  near-constant  volume  process  as 
long  as  the  energy  addition  is  limited.  When  the  energy  added  reaches  a  specific  critical  level, 
inertial  confinement  fails  and  the  entire  reaction  event  is  described  by  fully  compressible 
equations.  This  result  provides  quantitative  chemico-physical  conditions  for  the  validity  of  the 
Zeldovich  spontaneous  propagation  model. 

2.  Modeling  Concepts 

A  finite  spheroidal  volume  V"  of  spatially  uniform  reactive  gas  mixture  ( p'() ,  T0',  p'0)  with  an 
average  radius  R\  is  heated  by  absorbing  thermal  power  /y(J/kg.)  from  an  external  source  on  a 
heating  scale  t'H .  How  can  the  modeler  quantify  the  response  of  the  gas  volume  (the  near-field) 
to  the  energy  deposition,  E'  =  P't'u  (J/kg.)  ?  What  are  the  magnitudes  of  the  induced  motion  and 
thermodynamic  variable  changes?  How  does  the  near-field  response  affect  the  unheated  gaseous 
environment  (the  far-field)?  Thermo-mechanical  concepts  can  be  employed  to  address  these 
questions.  Kassoy  [2]  has  used  asymptotic  methodologies  to  describe  the  consequences  of 
energy  deposition,  E' »  e'Q  =  C’T^  (the  initial  internal  energy)  into  a  finite  volume  of  inert  gas 
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when  s  =  t'H  /t'A  «  1,  where  t'A  =  R'/ a'Q  is  the  characteristic  acoustic  time  and  a0  is  the  speed  of 

sound  based  on  the  initial  spot  temperature  7]'.  During  this  brief  period  of  heating,  an  acoustic 

disturbance  can  propagate  only  a  small  distance  compared  to  the  spot  dimension  R' .  When  E’  is 
less  than  a  specific  critical  value  the  heating  process  is  nearly  constant  volume  (inertial 
confinement),  characterized  by  asymptotically  small  values  of  density  variation  along  with  a 
large,  lowest  order,  constant  volume  temperature  increase  as  well  as  a  concomitant  large  rise  in 
pressure.  During  this  period  of  near-inertial  confinement,  the  thermo-mechanically  induced 

near-field  gas  Mach  number  (based  on  the  local  hot  speed  of  sound,  a'  =  sjyR'T'  is  very  small, 

implying  that  the  preponderance  of  the  absorbed  thermal  power  is  used  to  enhance  the  internal 
energy  of  the  gas.  Very  little  thermal  power  is  converted  to  kinetic  energy.  A  limited  range  of 
^'-values  exist  for  which  the  Mach  number  of  the  gas  expelled  from  the  edge  of  the  volume  is 
asymptotically  small.  The  “piston”  effect  of  the  expansion  generates  only  small  amplitude 
acoustic  waves  in  the  far-field.  Larger  values  of  £',  but  less  than  the  aforementioned  critical 
value,  are  associated  with  stronger  “piston”  effects  that  are  the  source  of  more  robust 
gasdynamic  disturbances  in  the  far-field,  similar  that  described  in  [22]. 

When  E'  reaches  the  critical  magnitude,  the  heat  addition  process  is  determined  to  be 
completely  compressible,  characterized  by  large  increases  in  temperature  and  pressure  as  well  as 
an  (2(1)  or  larger  internal  Mach  number.  There  is  significant  conversion  of  thermal  to  kinetic 
energy  in  the  near-field.  The  relatively  violent  expansion  process  is  the  source  of  strong  blast 
waves  in  the  far-field  [23]. 

The  application  of  these  thermomechanical  concepts  to  a  reactive  gas  is  affected  by  the 
highly  nonlinear  temperature  dependence  of  exothermic  chemical  reactions,  particularly  for  high 
activation  energy  kinetics.  As  a  result  there  will  be  consequential  differences  in  the  outcomes 
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arising  from  rapid  chemical  heat  release.  The  questions  raised  in  the  first  paragraph  of  this 
section  remain  germane  for  the  reactive  case.  The  answers  may  help  to  understand  how 
transient,  spatially  distributed  heat  release  in  a  reactive  flow  can  be  the  source  of  acoustic  and 
gas  dynamic  disturbances  (e.g.  solid  and  liquid  rocket  engine  dynamics,  internal  combustion 
engines). 

The  conceptual  model  description  provides  physically-based  insights  into  the 
thermomechanics  of  gases  affected  by  localized,  relatively  rapid,  transient  energy  addition.  The 
mathematical  model  used  in  Ref.  2  provides  quantitative  assessments  of  the  physics,  a  systematic 
scheme  based  on  formal  asymptotic  methods,  including  a  complete  understanding  of  the  model 
limitations.  A  well-defined  thermomechanical  analysis  for  a  reacting  flow  is  described  in  the 
next  section. 

2.1  Mathematical  Model 

Imagine  thermal  power  deposition  into  a  spheroidal  volume  of  characteristic  dimension  R' 
during  a  heating  time  period  t'H  such  that  t'COLL  t'H  <sc  t'A  -  R'/a'0  where  t'COLL  is  the  collision 
time  for  the  reactive  gas  mixture.  The  gas  is  initially  at  rest  with  density,  pressure,  and 
temperature/?',  P'}  and  7)',  respectively,  while  the  ambient  speed  of  sound,  a'Q  =  (yR'T^f^ .  A 

simple  one-step  high  activation  Arrhenius  exothermic  reaction,  R  — >  P,  is  assumed  to  describe 
the  chemistry. 

The  mathematical  model  is  based  on  the  equations  for  a  compressible  perfect  reactive  gas 
with  transport  effects  included.  Conservation  and  state  equations  in  three  spatial  dimensions 
(/,  0,  ^)and  time  (/')  can  be  written  in  the  nondimensional  form: 
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pT  +  eV  ■  pu  =  0 


(1) 


^  n 

p\_uT  +su  •  VmJ  =  -s - h  J3V 

7 


P 


(r~  0 


[Tt  +  su  •  VT  J  =  -£pV  ■  u  +  y AHepQ  +  /? 


r  c 

(y-l)  Pr 


+  yO 


(2) 

(3) 


Yt+£u-VY  =  -SYe t*'t+P--D  ,  p  =  pT  (4a,  b) 

Sc  p 

where  the  three-dimensional  vector  operators  are  in  spherical  coordinates. 

The  equations  can  be  nondimensionalized  initially  with  variables  (unprimed)  defined  by 


(p,p,T)  =  (p'/p'0,  p'/p'Q,  T'/tf)  ,  u  =  m'/ a'Q 

(5a, b) 

x-x'/R' 

T=t'/t'H 

(6a, b) 

ii 

^i 

C  =  C'/  k'°T°  ®  =  ®  /  //(,a<)  D  =  D'l ^ 

’  /  R’2  ’  /  R'2  ’  /  R'2 

(7a,b,c,d) 

Q=a(r,ff,f>,T)+nYeT‘/T  ,  TA=TJX  ,  T^E'/R' 

(8a,b,c) 

S  ft  A  ’ 

/S=4Av  .  s=,'hb' 

(9a,b,c) 

n=A H‘"S 

AH, 

.  •  ah,=ah;  14 

(10a,b,c) 

where  V,  C,  ®  and  D  are  symbols  representing  the  dimensionless  viscous  vector,  scalar 
conduction,  scalar  viscous  dissipation  and  scalar  mass  diffusion,  respectively.  In  addition, 
,  k'Q  and  22^  are  the  characteristic  viscosity,  conductivity  and  diffusivity  at  T'r  respectively. 

Equation  8a  contains  the  external  power  source  Qe  -  Q'J (AH'e/t'H )  where  Q'e  is  the  dimensional 
power  source  and  A H'clt'H  is  the  characteristic  dimensional  power  deposition,  composed  of  the 
characteristic  energy  A H'e  added  on  the  heating  time  scale  t'n.  In  addition,  Q  is  a 
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nondimensional  pre-exponential  factor  defined  in  (10a),  where  A Hch  is  the  nondimensional 
chemical  heat  of  reaction,  measured  with  respect  to  a'Q2  -  /(/-l)^,  a  surrogate  for  the  initial 
internal  energy,  e'Q  at  ro'  and  8  is  defined  in  (9c)  as  the  ratio  of  the  heating  time  to  the  inverse 
pre-exponential  factor  B'  \  for  a  one-step  reaction.  Similarly,  A He,  in  (10c),  is  the 
nondimensional  externally  added  energy,  also  measured  with  respect  to  a'2.  Y  is  the  mass 
concentration  of  fuel,  initially  distributed  uniformly  at  1.  The  nondimensional  activation 
temperature  Ta  is  defined  in  (8b)  with  respect  to  the  initial  temperature  T’v  where  the 
dimensional  activation  temperature  is  defined  in  (8c)  in  terms  of  the  dimensional  activation 
energy  E'{)  and  the  universal  gas  constant  R' .  It  is  assumed  that  the  characteristic  heating  time 

t'H  <sc  t'A  —  R'/a'0  so  that  the  parameter  s  defined  in  (9a)  is  very  small.  It  is  noted  from  (9b)  that 
the  characteristic  viscous  time  t'v=R'2/v'Q  is  much  larger  than  the  heating  time  because 


p  =  t'H/t'v  =  £7  0 1  where  the  acoustic  Reynolds  number  a'QR'/v'Q  is  large  for  R' :»  O(l0  8m). 


The  Prandtl  and  Schmidt  numbers  as  well  as  y  have  standard  definitions  and  are  treated  as  0(1) 
constants,  relatively  to  the  small  parameter  s. 

Equation  (3)  is  written  for  a  constant  specific  heat  to  provide  a  relatively  elementary  model. 
This  assumption  has  the  disadvantage  of  precluding  activation  of  internal  energy  modes  at  high 
temperature.  Limitations  are  discussed  in  Section  3. 

The  rest-state  initial  conditions  are 


r  =  0  :  p  =  P  =  T  =  1  ,  u  =  0  ,  Qe  =  0 


(11) 


In  general,  the  spatially  distributed,  time-dependent  external  source  function  can  be  described  by 
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(12) 


Qe{r,  0,  <p,  t)>0  ,  r<r.(0,  <p) 

Qe(r ,  0,  <p,  t)  =  0  ,  r>rs(9 ,  <p) 

where  r  (Ojp'j  =  r’(0,(p'j/R'  describes  coordinates  of  the  outer  surface  of  the  heated  region.  The 

volume  beyond  the  spot  surface  is  not  heated  by  the  source. 

The  nondimensional  equations,  (l)-(4),  describe  reactive  flow  dynamics  for  a  wide  range  of 
length  and  time  scales.  Table  1  provides  estimates  for  the  acoustic  t'A  and  viscous  t[  times  for 
several  values  of  R'  and  the  allowable  range  of  heating  time  values  t'H,  such  that 
t'coLL  ^  r',\  ■  The  numerical  values  for  t'A  are  defined  for  a  gas  initially  at  T’{]  -  300  K  and 
p'a  =  1  atm.  These  ambient  values  are  associated  with  a  mean  free  path  A’  -  5.3 x  10  s  m,  a 
collision  time  t'COLL  =  1 ,6x  10  10  s  and  a  characteristic  kinematic  viscosity  v'  =  1 .64  x  10  5  m2/s. 

Asymptotic  methods,  based  on  the  time-ratio  limit  s  — >  0,  are  used  to  obtain  reduced 
equation  systems  for  several  phases  of  the  heating  process,  reaction  initiation,  and  the 
accompanying  generation  of  gasdynamics  waves.  It  is  assumed  that  Pr  =  0(1),  Sc  =  0(1),  and 
that  (3  — >  0  (see  7b)  because  the  viscous  diffusion  time  defined  previously  is  typically  orders  of 
magnitude  larger  than  the  heating  time.  Finally,  it  is  assumed  that  the  activation  temperature 
defined  in  (8b, c),  TA  -  TA  (y;) :»  1,  where  the  specific  functional  dependence  is  defined  later  in 
the  analysis. 
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Table  1  -  Characteristic  values  of  the  acoustic  time  t'A  and  the  viscous  time  t'v  and  the  allowable 
range  of  heating  times  t'H  for  a  range  of  characteristic  spot  dimensions  R'  when  T’}  -  300  K. 
Pq  - 1  atm  and  v'0  =  1 .64x1 0  5  m2s 


R'(  m) 

t'A  =  R'/a'0(s) 

4(s) 

t'v=R'2/v'0(  s) 

10-2 

3  x  10'5 

<10'6 

6.1 

10"3 

3  x  10~6 

<10~7 

6.1  x  10'2 

10"4 

3  x  10'7 

<10'8 

6.1  x  10'4 

10'5 

3  x  10'8 

<10'9 

6.1  x  10'6 

10'6 

3  x  10~9 

tf 

^COLL 

6.1  x  10'8 

2' 

tf 

lCOLL 

— 

- 

2.2  Inert  Heating  Period:  t'  =  0(t'H  ) 

(a)  Internal  spot  dynamics:  r  <  rs  ( 6 1  (p) 

The  limiting  forms  of  (1),  (2)  and  (4a)  for  a  — >  0  and  TA  (s)  — >  qo,  pT  =  u  T  =  Yz  =  0  and  the 
relevant  initial  conditions  in  (22)  imply  that  uq  =  0  and  /9o  =  yo  =  1,  where  the  zero  subscript 
represents  the  lowest  order  asymptotic  approximation  to  each  of  the  dependent  variables.  It 
follows  from  (4b)  that  po  =  To  where  the  temperature  distribution  is  described  by  the  reduced 
version  of  (3)  and  the  appropriate  initial  condition  respectively; 
dT 

— ~  ( r,0,<p,r ) 

dr  V  ;  ’  (13a, b) 

r  =  0  ;  T0  =1 

Equation  (13)  describes  constant  volume  heating  due  to  energy  deposition  represented  by  Qe. 
The  characteristic  changes  in  temperature,  AT'/T^  -  O(l)  and  pressure,  A p'Q/  p'Q  =  O  ( 1 ) ,  occur 
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for  A He  =  0(1)  and  finite  values  of  z.  Chemical  heating  plays  no  role  in  this  time-frame  because 
the  activation  temperature,  TA  in  (8b, c),  is  assumed  to  be  large. 

The  fluid  dynamic  response  to  the  spatially  distributed  temperature  and  pressure  (To,  po)  can 
be  obtained  by  substituting  the  transformed  density  and  velocity  variables, 

p  =  l  +  s2R  ,  u  =  sU  (14a,b) 

Into  (1),  (2)  and  (11)  and  taking  the  limit  £— >  0.  The  reduced  momentum  and  conservation  of 
mass  equations,  along  with  initial  conditions, 


U0r 

© 

i 

ii 

,  Uo(r,0,q>,  0)  =  0 

(15) 

r 

=  -v-t/0 

,  RQ(r,6,(p,0)  =  0 

(16) 

describes  gas  motion  induced  by  the  pressure  (temperature)  gradient  created  by  spatially 
distributed  heating  from  Qe,  and  the  concomitant  density  variation.  Equations  (5a)  and  (14b) 

can  be  used  to  show  that  the  local  gas  speed  is  characterized  by  (9  (fizz')  or  the  local  Mach 

number,  ulT  “  =  O(s),  is  very  small.  Solutions  to  (13),  (15)  and  (16)  can  be  written  in  quadrature 
form; 


To=l  +  r(r-l)^He\Qe(r,0,<P,f)dT 

0 

II 

o 

a. 

(17a, b) 

1  r 

U,  =  —  [V7Itfr 

(18) 

rl 

T 

R0  =  -  {  V  •  U0d  x  , 

(19) 

o 

for  r<rs(6,  (p)  defined  in  (12)  where  V  represents  the  gradient  operator  in  spherical  coordinates. 
Relatively  elementary  (analytical)  solutions  can  be  obtained  for  a  model  source  function, 
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Qe=w(r)F(0,q>)g(  r)  , 

(20a) 

where 

w(r)>0  ,  r<rs{9,(p )  , 

(20b) 

w(r)  =  0  ,  r>rs(0,(/))  , 

(20c) 

o 

ll 

o 

&0 

(20d) 

g(r)>0  ,  r  >0  , 

(20e) 

limg(r)  =  goo  >0  ,  g^  =  constant 

r-»oo 

(20f) 

The  energy  deposition  associated  with  the  source  is  confined  to  a  targeted  region  r  <  rs  ( O.cp),  and 
is  assumed  to  asymptote  to  a  purely  spatial  distribution  for  large  values  of  the  time  variable  r. 

The  time-varying  spatial  variation  in  temperature,  found  by  substituting  (20a)  into  (17),  is 
determined  by  the  product  w{r)F{0,(p)\ 

T 

To(r,0,(p,T)  =  l  +  y(y-l)AHew(r)F(6,<p)^  g(f)dt  (21) 

0 

Then  from  (18)  and  (19)  the  velocity  and  the  density  perturbation  are  described  by 

T  T 

U0{r,6,(p,r)  =  -{y -l)  AHeV[w(r)F  (6>,^)]  JdfJ  g(a)d(T  (22) 

0  0 

T  T  <7 

RQ{r,d,(p,v)  =  (y-l)AHiy2\_w(r)F(d,(p)]^dt^d(j^  g((j)d(j  (23) 

0  0  0 

The  definition  of  the  source  function  w(r)  in  (20b  ,c)  can  be  used  in  (21)  and  (22)  to  show  that  at 
the  edge  of  the  spot  r  =  rs  (O.cp),  the  temperature  To(rs,d,(p,  r)  =  1,  the  radial  velocity  component  is 

Uor(rs’0’<P’T)  =  -(r-1)AHeF(0,<p)— Vj^-\dT^g{(j)d(j  (24) 

ar  o  o 
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while  the  angular  components  Uoe  and  Uo,p  vanish.  Equation  (23)  shows  that  the  density 
perturbation  at  the  edge  Ro(rs, 0,  (p,  z)  can  be  nonzero.  When  dw/dr(rs)  ^  0,  the  radial  velocity  at 
the  spot  edge  in  (24)  is  the  source  of  acoustic  waves  in  the  adjacent  cold  gas  environment  that  are 
described  in  section  2.3,  below. 

Example  solutions  of  (21)-(23)  for  a  simple  spherically  symmetric  source  defined  by 
w(r)  =  cos(nr/2)  ,  F(G,(p)  =  1  ,  g(r)  =  g„(l-^r)  ,  r<  1  (25a, b,c) 

are 


T0  =l  +  y(r-l)Miegxco^[z-l+e-T) 


(21a) 


U0r  =(y-l)AHegJ^sm^- 


V  A 

- — z-e~T  +  1 

2  j 


(22a) 


R0=-(Y-l)AHeg^ 


n  nr  .  nr 
—  cos - b  srn  — 

v  2  2  2  j 


V  r2  ,  . 

- br-l  +  e 


(23a) 


where  the  spatial  structure  is  determined  by  w(r )  in  (25a)  and  the  time  dependence  follows  from 
g(r)  in  (25c).  The  temperature  distribution  in  (21a)  is  radially  symmetric  with  a  maximum  at 
each  value  of  r  located  at  r  =  0,  and  a  minimum  at  r  =  1.  The  spatial  amplitude  at  each  radial 
location  grows  linearly  with  r»  0(1)  due  to  the  sustained  energy  input.  Meanwhile,  the 
thermomechanically-induced  radial  fluid  speed  increases  from  zero  at  r  =  0  to  a  maximum  at 


r  =  1 ; 


r  =  l  ;  U0r  = 


C  n^ 


(r~ !) 


1 

- r+l-e 

v  2  , 


(26) 


When  t'  =  0(t'H )  or  r=  (9(1),  this  accelerating  edge  speed,  associated  with  an  0(e)  local 


Mach  number,  given  the  definition  in  (14b),  is  a  “piston-like”  source  of  acoustic  waves  [3].  The 
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waves  can  move  a  distance  ArJ,  =  O(t'Ha'0 )  =  O(eR'),  small  compared  to  the  spot  dimension 
itself.  It  should  also  be  noted  that  the  fluid  itself  is  convected  a  characteristic  distance 
A  r'F=o(e2R'). 

The  mass  loss  at  the  edge  (r  =  1)  compatible  with  (26)  causes  a  density  decrement  within  the 
hot  spot,  described  by  (23a).  The  spatial  structure  function  is  positive  for  O  <  r  <  1  and  the 
spatial  amplitude  increases  like  0(  z3)  for  r  »  0(1).  Note  that  Ro  (r  =  1,  r)  <  0  implying  that  the 
density  (14a)  is  slightly  smaller  than  the  initial  value. 

2.3  Local  Linear  Acoustics:  r  >  1 

The  acoustic  region  variables,  defined  by 

r-l  +  er  ,  r> 0  ,  u=eU  ,  (p— 1,  P~h  T -l)  =  e(p,R,T^  ,  (27a,b,c) 


with  Q  -  0  are  used  to  describe  an  acoustic  disturbance  in  a  thin  boundary  layer  around  the  edge 
of  the  spherical  spot  for  the  example  source  functions  in  (25).  The  0(e)  magnitude  of  the 
thermodynamic  variables  perturbation  is  compatible  with  an  0(e)  edge  Mach  number  defined  by 
(26)  and  (27b),  valid  when  r  <?C  0( l/e112).  Equation  (27)  can  be  used  in  (l)-(4)  and  (11)  to  find 
linear  acoustic  equations  valid  in  the  limit  e  — >  0, 


R0  +V  -U0  =0  ,  U=—VpQ 

r 


(28a,b) 


A,  =  rK  ,  Po=To+Ro 


(29a,b) 


Equations  (28)  and  (29)  can  be  combined  to  find  a  simple  wave  equation  for  the  radial  velocity 
component  Uor 

d2Un 


dr 


=  V  U 

2  V  U  Or 


(30) 
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where  (p  is  the  function  on  the  right  side  of  (26)  with  r  replaced  by  1  -  //.  The  acoustic  velocity 
disturbance  described  by  (32),  and  the  associated  thermodynamic  variations  found  from  (28a) 
and  (29)  are  smooth  functions  with  benign  front  properties  but  are  valid  only  for  77 0(\l£V1). 
The  limitation  arises  because  the  boundary  condition  in  (26)  defines  an  accelerating  “piston” 
speed  that  eventually  exceeds  the  small  Mach  number  condition  implicit  in  the  derivation  of  (28) 
and  (29).  Uqt  and  its  first  spatial  derivative  are  zero  at  the  front  77  =  1.  The  discontinuity  with 


respect  to  the  undisturbed  value  for  77  >  1  is  in  the  second  spatial  derivative.  These  function 
characteristics  arise  from  the  properties  of  the  time-resolved  source  in  (25c).  One  may  note  from 
(26)  that  the  effective  “piston  speed”  associated  with  the  expanding  hot  gas  is  O(t')  for  r<tcl. 
This  relatively  gentle  startup  process  should  be  compared  with  typical  instantaneous  acceleration 
problems,  which  display  solutions  with  discontinuous  functions  at  the  front  [22].  In  contrast,  the 
edge  speed  increases  like  0{t)  for  large  time  values,  implying  that  successively  stronger 
disturbances  will  appear  in  the  far-field  as  time  evolves.  It  is  noted  from  (27c)  that  the  acoustic 
temperature  change  is  0{s)  implying  that  the  chemical  reaction  is  suppressed  in  the  external 
environment  for  r=  (9(1). 
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Finally,  it  is  noted  that  more  general  acoustic  results  can  be  obtained  by  using  the  edge 
speed  in  (24)  as  the  “piston-like”  source  of  disturbances.  The  acoustic  characteristics  will 
depend  on  the  spatial  and  time  dependence  of  the  model  source  function  in  (20a). 

2.4  Transport  Effects  near  r  =  1. 

The  rising  temperature  in  the  heated  region  caused  by  the  source  defined  in  (12)  will,  in 
general,  lead  to  a  finite  temperature  gradient  near  the  edge  of  the  spot  and  conduction  into  the 
cold  environment.  An  assessment  of  the  conduction  process  on  the  heating  time  scale 
t'H(r  =  O(l))  can  be  obtained  by  transforming  the  energy  equation  (3)  with  scaled  variables  that 

describe  small  temperature  variations  from  the  edge  value  T  =  1  in  a  thin  boundary  layer  in  the 
vicinity  of  rv(  0,  (p).  A  balance  of  energy  accumulation,  the  absorbed  heat  and  the  conduction  term 
can  be  used  to  show  that  T  -  1  =  0(f:lU2)  and  r  -  rs(0.(p)  =  0( /i1/2).  Given  the  definition  of  /?  and  s 
in  (9a, b),  it  follows  that  ft  can  be  less  than  e.  In  the  asymptotic  limit,  the  conduction  boundary 
layer  is  much  thinner  than  the  acoustic  layer  defined  in  (27a)  and  the  local  temperature  variation 
0(J3  )  is  minute  relative  to  the  O(s)  acoustic  temperature  disturbance.  Finally,  the  ratio  of  the 
characteristic  heat  loss  rate  at  the  surface  of  the  spot  to  the  volumetric  power  absorption  is  0{/3), 
which  is  negligible  in  the  limit  s  — >  0. 

2.5  Asymptotic  Solution  Properties:  r— »oo 

A  study  of  the  solution  properties  for  z  — »  oo  from  (17)-(19)  provides  an  opportunity  to 
describe  how  the  rising  spot  temperature  can  initiate  a  high  activation  energy  reaction  associated 
with  the  kinetics  term  e~TA'T  where  limT^  (f)  — »oo.  When  the  exponential  term  is  (9(1),  or  T  = 

0(Ta),  one  can  anticipate  a  strong  exothermic  reaction. 
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Asymptotic  behavior  of  the  temperature,  pressure,  speed  and  density  can  be  obtained  easily 
for  the  generalized  deposition  heat  source 

Qe  =  f(r>0><P)8(T)  (33a) 

where 


f(r,0,(p)>  0 

,  r<r(0,(p) 

(33b) 

f(r,0,(p)  =  0 

,  r>rs(e,cp) 

(33c) 

and  g(z)  has  the  properties  defined  in  (20e)  and  (20f).3 

Equations  (17)-(19)  can  be  evaluated  to 

show  that 

limro  ~  r  , 

r— >00 

lim  pQ  ~  r 

r— »  oo 

(34a, b) 

limf/0  ~  r2 

r— »oo 

lim  Rlt  ~  r3 

r^oo 

(34c, d) 

In  addition,  the  local  Mach  number 

in 

ii 

3/2 

1  ~  ST 

(35) 

The  high  activation  energy  reaction  is  initiated  when  To  ~  r=  0(73).  Then  (14),  (34),  and  (35) 
imply  that 

T  =  0  (Ta  )  ,  p  =  0(TA )  ,  u=0(sT2a)  (36a, b,c) 

p-\  =  0(s2Tl)  ,  M  t  =  o{sTa2)  (36d,e) 

where  TA  »  0(1).  The  dependence  of  the  local  Mach  number,  Me,  on  both  s  and  Ta  can  be  used 

to  discriminate  between  two  kinds  of  evolutionary  reactive  gas-dynamic  processes  on  an 
extended  time  scale: 


3  Other  forms  of  time-dependence  can  be  evaluated. 
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(a)  incompressible  limit:  Mf  «:  1  ,  1  <scTA  «;  2/3 )  (37a, b) 

p  =  T  <^o(s  213)  ,  m«o(^"1/3)  ,  (p-\)  =  o(\)  (37c,  d,e) 

(b)  compressible  limit:  Me  =  <9(l)  ,  TA  =o( c  2/3  j  (38a,b) 

p  =  T  =  o(£-2'3)  ,  m=o(^1/3)  ,  (p-l)  =  0(l)  (38c, d,e) 


Equations  (37b)  and  (38b)  provide  explicit  limitations  on  T\  for  a  prescribed  s,  or  vice  versa.  The 
incompressible  case,  characterized  by  a  small  local  Mach  number  (37a)  and  small  density 
variation  (37e),  can  be  expected  to  retain  the  near  inertial  confinement  properties  found  for  the 
initial  heating  period  r  =  0(  1).  In  contrast,  the  compressible  limit  properties  in  (38)  imply  a 
more  significant  transient  gas  dynamic  process,  including  the  conversion  of  significant  amounts 
of  thermal  energy  to  kinetic  energy. 

A  physical  explanation  for  the  failure  of  inertial  confinement  when  Ta  =  0{s '  )  can  be 
obtained  by  considering  the  ratio  of  the  dimensional  extended  heating  time  t’He  —t'HTA ,  derived 

from  r=  0(TA),  to  the  dimensional  local  acoustic  time  t'Ae  =  R'/a'  =  t’Al\fr  =  o{t'A!TA2\ 


^ He  „y3/2 

t'  A 

lAl 


(39) 


Inertial  confinement  can  prevail  only  when  the  ratio  is  small,  corresponding  to  case  (a).  When 
Ta  =  0(  A2/a),  the  ratio  is  0(1),  implying  that  acoustic  wave  can  traverse  the  spot  on  a  time  scale 
comparable  to  that  for  extended  heating. 

The  scaling  contained  in  (36)-(38)  can  now  be  employed  to  define  rescaled  variables  and 
new  limiting  equations  for  the  reactive  initiation  portion  of  the  process. 
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3.  Reaction  Evolution:  t'  =  0(TAt'H) 

3.1  Incompressible  Limit 

The  temporal  nonuniformities  of  the  inert  heating  period  solutions  in  (34)-(36)  imply  that  the 
rescaled  variables  for  the  reactive  incompressible  limit  are; 


t  =  Tat  ,  l«rA  «o(^-2/3)  (40a,b) 

(p,  T)  =  (tap,  TaT )  (40c, d) 

p  =  l  +  s2T3R  (40e) 

u  =  sT2u  (40f) 


These  variables  describe  an  extended  period  of  thermal  energy  deposition  from  a 
combination  of  external  source  and  exothermic  chemical  heat  release.  The  definition  of  Qe 
below  (5)-(10)  implies  that  the  characteristic  dimensional  power  deposition  during  the  inert 
period  r=  (9(1)  is  given  by  Q'e  =0(AH'e/t'H)  =  0(a'2/t'H \  =  O(e'0/t'H)  where  definitions  below 

(5)-(10)  have  been  used.  It  follows  that  energy  deposition  during  that  period  t'  -  0(t'H )  is 
simply  Q'et'H  -  0(e’>‘).  In  contrast,  energy  deposition  on  the  extended  time  scale  t'  -  0(TAt'H  )  is 
Q'et'HTA  =  O ( ef0TA  )-0( C'T'a  ) .  The  inequality  in  (40b)  implies  that  TA  <§:  ro' / s213 .  It  follows  that 
characteristic  energy  deposition  for  the  incompressible  limit  is  restricted  to 
Q'et'HTA  «  o{cXl£2,3)  =  0(e'ol£213)-  I11  other  words,  the  nearly  inertially  confined  heating 

process  will  occur  only  if  the  energy  added  during  the  extended  heating  period  is  sufficiently 
small.  This  class  of  limitation,  identical  to  that  found  in  Ref.  [2]  for  inert  gas  thermomechanics 
must  be  a  universal  characteristic  of  rapid  localized  heat  addition  to  a  variable  density  fluid.  The 
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extreme  case  of  energy  addition  corresponding  to  TA  =  0(  \hrB)  is  associated  with  a  fully 
compressible,  high  internal  Mach  number  flow.  Details  are  given  in  Section  3.2. 

The  temperature  and  pressure  changes  are  large,  while  the  density  variation  is  small. 
Examination  of  (40f)  shows  that  the  nondimensional  speed,  based  on  a'0  (see  (5b)),  can  take  on  a 

wide  range  of  values.  In  particular,  at  the  edge  of  the  spot  (T"  -  V),  the  edge  Mach  number 


«'/«;«  o(i) 

,  1  «  Ta  «  £  1/2 

r«e  1/2 

(41a) 

K/a'0=o(  f) 

II 

ro 

t  =  0(s-1'2) 

(41b) 

«!K»o(i) 

,  £■ 1/2 « 7; « £■ 2/3  , 

,  T  <SC  S~2B 

(41c) 

where  u'e  =  I //I  and  the  third  inequality  in  each  of  (41)  follows  from  the  time  transformation  in 
(40a).  The  subsonic  edge  Mach  number  in  (41a)  is  compatible  with  the  linear  acoustic  solution 
in  (32a),  due  to  the  specified  time  extent  r  <sc  £  l/2. 


Accordingly,  the  character  of  the  gasdynamics  in  the  cold  external  environment,  driven  by 
gas  expansion  at  the  edge  will  differ  considerably,  depending  on  the  particular  value  of  TA.  In 
contrast,  the  local  Mach  number  inside  the  spot  Me  =  \j/\/a'  =  \ij\j *Jt  =  < tc  1. 

Equation  (40)  can  be  substituted  into  (l)-(4)  to  derive  the  describing  equations  for  the 
extended  heating  period  r=  0(TA). 


Rf  +V \pu)  =  0 


p  (Q,  +  e-Tl  (a  •  v)  a )  =  -  ^  -  pr^v 


(42a) 

(42b) 


P  (ff  +  e-Tlu  •  Vf )  =  sXpV  •  u  +  yAHePQ  +  J3T1A 


7  °  +y£2Tfo 


(y  —  l)  Pr 


(42c) 
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1  IT  ,  P  — •  D 


Y.  +  slT*u  •  VY  =  -8TAYe  +  T”  — 

Sc  p 


(42d) 


Q  =  Qe+DYe 


(42e) 


P  =  [\  +  £2TIr)t 


(42f) 


where  power  law  dependence  for  viscosity,  conductivity  and  diffusion  coefficient,  respectively, 
defined  by  ju  =  T\  k  =  T'  and  D  =  Tnlp,  where  n  <  1  and  m  <  1,  are  assumed  for  the 
nondimensional  transport  equations  to  provide  a  model  that  includes  high  temperature  transport 
effects. 

Eq.  (42)  reduces  to  a  relatively  elementary  system  in  the  limit  s  — >  0  because  all  of  the 
transport  term  parameters  are  vanishingly  small  (see  the  definition  of  (3  in  (9b),  and  Table  1  for 
characteristic  values  of  t'v,  while  TA  <50  is  typical). 


R0l  +  ^  ’  Md  0 


VA> 

r 


(43a,b) 


(r- !) 


T0f  =  yAHe  2e+ny//f»  ,  Qe  =  Qe  ( r,0,<p,f ) 


(43c, d) 


4  =-STAY()e-uf<>  ,  p0=T0  (43e,f) 

where  Q  and  S  are  defined  in  (10a)  and  (9c),  respectively,  and  Qr  represents  the  effect  of 

continued  external  thermal  power  deposition.  It  is  assumed  that  the  heat  of  reaction  A Hch  = 
0(Ta )  in  order  to  assure  that  the  temperature  increase  due  to  exothermicity  is  of  the  same 
magnitude  as  that  due  to  external  source  effects  on  the  time  scale  f  =  0(1).  In  addition,  it  is 
assumed  that  8  =  0(l/TA )  to  assure  0(1)  reactant  consumption  on  the  time  scale  f  =  0(1).  Both 

assumptions  are  physically  viable.  Equation  (43)  describes  an  extended  period  of  nearly- 
inertially  confined  thermal  power  deposition,  including  a  small  change  in  density  in  an  internally 
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low  Mach  number  expansion  process.  Equations  (43c)  and  (43e)  can  be  combined  to  define  a 
Schwab-Zeldovich  variable  [24]  used  to  find  F0  =  Y0(f0); 


Yn=\  +  ^U-]Qedf 


AH 


ch  0 


^chy{r-]) 


ah„=^A=o(  i) 


(44) 


Equation  (44)  can  be  used  to  find  the  maximum  temperature  Tq(t  =  oo)  associated  with  complete 


reactant  consumption  Yq  =  0 


T0(f=co)  =  r(r-i) 


AH 


ch 


■A He^Qe(r,6,(p,  f)df 


(45) 


where  the  first  term  on  the  right-hand  side  represents  the  adiabatic  temperature  rise  from  constant 
volume  complete  reactant  consumption  and  the  second  term  is  the  total  temperature  rise 
associated  with  the  sustained  external  source  energy  deposition. 

An  equation  for  f  is  found  by  combining  (43c),  (44)  and  (45) 


f0,=r(r-i)w,Q,+-^ne-“f- 

AHch 


(46) 


The  initial  condition  on  f0  is  found  formally  by  using  the  matching  expression  from  (40d) 


T0(f  — >0)  =  T(j  — >oo)/rA  =0,  meaning  that  looking  back  in  time  to  the  period  r  =  (2(1),  the 
temperature  is  asymptotically  small  compared  to  that  when  f  =  0(  1).  Equation  (46)  can  be 
integrated  numerically  to  find  To(r,0,(p,f )  for  a  specified  external  heat  source  Qe.  Reactant 

consumption  is  defined  by  (44).  The  locally  high  pressure,  defined  by  (43f),  can  be  used  in  (43b) 
to  find  the  local  thermomechanic  ally  induced  fluid  motion.  Finally  the  density  distribution  is 
found  by  integrating  (43a).  The  equation  system  bears  a  striking  resemblance  to  that  discussed 
in  Ref.  [2]  in  the  context  of  the  response  of  an  inert  gas  to  extremely  rapid  transient,  spatially 
resolved  energy  addition. 
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The  systematically  derived  reduced  equations  (43-46)  provide  a  formal  validation  for  the 
Zeldovich  [19]  detonation  initiation  model  based  on  constant  volume  thermal  explosion 
evolution.  The  current  analysis  provides  specific  parametric  conditions  for  which  the  model  is 
valid  and  those  for  which  it  is  not  (see  Section  ?  ).  Related  concepts  have  been  developed 

earlier  by  Jackson  et  al.  [9]  and  Short  [10,  11]. 

Once  f0  and  To  are  found  from  (46)  and  (44),  respectively,  43b  and  f  can  be  used  to  find  the 
thermomechanically-induced  velocity 


Finally,  the  small  density  distribution  represented  by  R0  in  the  context  of  (51e)  is  found  from 

f 

R0  =  -  j  V  •  u^d  f  (48) 

o 

A  numerical  solution  to  (46)  for  T0(r,f)  with  the  initial  condition  T0(r,0)  =0  has  been  found  for 
the  axisymmetric  source  Qe  =  (l- f)cos;r(r/2), (0  <  r  <  1)  when  0<f<l,  and  <2f=0when 
f  >  1  for  a  few  values  of  the  parameters  A He  and  Q  with  A //,/,=  1  (see  (44)).  Beyond  the  edge  of 
the  spot  located  at  r  =  1,  <2  =  0.  The  radial  dependence  of  the  scaled  temperature  f0  is  given  in 

Fig.  1  for  several  values  of  the  scaled  time  when  A He  =  5  and  Q  =  10.  At  each  value  of  r,  the 
temperature  varies  from  0  to  a  maximum  value  obtained  from  (45).  Source  heating  causes  the 
temperature  to  rise  in  0  <  r  <  1  for  f  <  0. 1.  Thereafter,  localized  chemical  heat  release  causes  a 
relatively  rapid  temperature  increase  in  the  hottest  portion  of  the  spot.  By  the  time  f  =  1  when 
the  source  is  turned  off,  approximately  85%  of  the  spot  has  reached  the  maximum  temperature 
value.  A  substantial  radial  temperature  gradient  appears  in  the  vicinity  of  the  r  =  0.9  and 
becomes  steeper  as  time  increases.  Fig.  2  shows  the  evolution  of  the  gradient  near  the  edge  of  the 
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Fig.  1.  Temperature  Profiles:  AH  =  5.  Q  =  10 


Fig.  2.  Temperature  Profiles:  AH  —  5.  Q  =  1 0 
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spot  for  much  larger  times.  The  temperature  distribution  is  evolving  to  a  discontinuous  function 
near  the  edge,  r  — »  1.  This  essentially  constant  volume  process  is  characterized  by  p0  =  T0,  so 

that  the  high  temperature  gas  region  is  also  at  a  locally  large  pressure.  Equation  47  implies  that 
the  large  gradient  region  is  a  source  of  localized  fluid  acceleration. 

The  narrow  spatial  region  denoted  by  the  large  gradient  separates  hot  burned  gas  on  the  left 
from  cold  unbumed  material  on  the  right.  Results  for  the  spatial  distribution  of  reactant  Y, 
shown  in  Fig.  3,  can  be  interpreted  in  terms  of  a  reaction  front  propagating  through  the  spot.  If 
one  examines  the  change  in  the  location  of  the  Y  =  0.5  value  with  time  during  the  period 
0.3  <f  <land  converts  the  results  to  dimensional  values  (e.g,  using  (6a,b),  (40a, d))  it  can  be 
shown  that  the  front  propagation  speed,  relative  to  the  local  (hot)  speed  of  sound,  a',  is 

Vp/a'  =  Vp/yjr/f  =  0(1)  with  respect  to  the  limit  s  0.  Results  in  Fig.  2  for  large  values  of  time 

shows  that  the  front  speed  decelerates  drastically  as  it  moves  into  the  colder  gas  ahead  where  the 
local  “explosion”  time  is  much  longer  than  that  in  the  hotter  gas  region.  Near  r  =  1“,  the 
temperature  distribution  is  essentially  frozen  for  f  >  1  because  the  source  is  extinguished  and  the 
gas  is  too  cold  to  produce  any  significant  amount  of  thermal  energy.  The  local  finite  gradient 
implies  that  fluid  expelled  from  the  hot  spot  edge  (e.g,  (47))  will  be  a  continuing  source  of 
acoustic  or  stronger  disturbances  (e.g,  see  (41))  in  the  external  cold  gas  with  consequences 
similar  to  these  described  in  [22]. 
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1.2 


Fig.  3.  Concentration  Profiles:  AH  =  5.  Q  =  10 

Temperature  results  for  the  parameter  values  A He  =  1  and  Q  =  50  are  given  in  Fig.  4.  In  this 
case  the  reaction  process  is  initiated  more  slowly  because  the  source  strength  is  reduced  by  a 
factor  five.  However  the  reactive  phase  is  stimulated  by  a  much  larger  value  of  Q.  The 
thermomechanic  ally  induced  radial  fluid  velocity,  shown  in  Fig.  5  for  a  limited  range  of  the  time 
variable,  is  caused  by  the  large  local  gradient  in  temperature  and  pressure  as  defined  by  (43b) 
and  (47).  Larger  time  values  are  displayed  in  Fig.  6.  These  results  show  rather  substantial 
growth  in  the  maxium  speed  values,  nondimensionalized  with  respect  to  the  cold  speed  of  sound, 
a'0.  The  maximum  values  listed  in  the  inset  in  Fig.  6  are  quite  large.  These  results  should  be 
interpreted  with  the  knowledge  that  the  local  Mach  number,  with  reference  to  the  hot  speed  of 
sound  is  very  small,  as  noted  below  (41c). 
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Fig.  4.  Temperature  Profiles:  AH  =  1 .  Q  =  50 


Fig.  5.  Velocity  Profiles:  AH  =  I .  Q  =  50 
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It  may  be  tempting  to  associate  the  pressure,  temperature  and  reactant  concentration  jumps 
observed  in  the  figures  with  a  reactive  shock.  However,  on  the  reactive  heating  time-scale  of 
interest  (40a)  the  physical  process  is  occurring  in  a  primarily  incompressible  fluid  and  the  front 
is  observed  to  stagnate.  In  fact,  an  ephemeral  hot,  high  pressure,  low  speed,  burned  out  spot  has 
been  generated  by  a  momentary  inertial  confinement  process  during  the  extended  chemical 
heating  time.  The  huge  pressure  gradient  across  the  front  will  be  the  source  of  a  localized 
relaxation  process,  characterized  by  hot  spot  gas  expansion  on  the  longer  acoustic  time  scale  of 
the  spot.  The  front  will  be  driven  into  the  surrounding  unburned  reactant  leading  to  the 
appearance  of  new  localized  combustion  waves.  Elucidation  of  this  relaxation  process  is  the 
subject  of  a  current  investigation. 

The  nondimensional  results  in  Figs.  1-6  can  be  interpreted  in  terms  of  a  dimensional 
example  based  on  the  following  data. 
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R'  =  103m  ,  ro'  =  300K  ,  p'=latm  ,  ^=340  m/s  ,  t'A  =2.9x10  6s 


t'H  =  10  7  s  ,  s  =  .034  , 

AHe  =  <9(105J/kg)  ,  TA  =  l/£112  =5.42  , 

where  the  TA  value  is  midrange  in  (40).  The  characteristic  high  temperature  TAT('t  is  about  1600K 
and  the  corresponding  high  pressure  TAp'Q  about  5.5  atm  (see  (40)).  The  small  density  change  is 
0(sm)  and  the  local  Mach  number  is  0(£1/4). 

3.1.1  Linear  Acoustic  Wave  Evolution:  t  -  (9(1) 

Acoustic  disturbances  of  O(s)  amplitude  described  in  Section  2.3  by  (27)  propagate  at  a 
characteristic  dimensional  speed  a'0  during  the  interval  t'  =  (9(7', ).  The  characteristic  distance 
traveled  is  O(eR'),  small  compared  to  the  hot  spot  dimension  R' .  During  the  extended  heating 
time  period,  defined  by  (40a),  the  wave  front  moves  a  greater  distance  characterized  by 
0{sTaR').  Meanwhile,  larger  linear  acoustic  disturbances  are  generated  by  the  small  edge  Mach 
number  defined  in  (41a)  during  the  time  period  characterized  by  r-Oi\).  The  variable 
transformations 


(p-\)  ,(p-\)  ,(T-\)  =  sT2A(p,pj)  ,  u  =  sT2au  (49a,b,c,d) 

r  =  TAf  ,  r  =  l  +  sTA^  (50a,b) 

and  the  limit  s  — >  0  can  be  employed  in  (l)-(4)  to  derive  linear  acoustic  equations  describing 
larger  acoustic  disturbances 


pf+V-u  =  0 


(51a) 
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(51b) 


(51c) 


(r- 0 


=  -V  •  u 


p  =  p  +  t  +0(s) 


(5  Id) 


where  the  differential  operators  are  written  in  spherically  symmetric  form 


These  equations  are  subject  to  initial  conditions  associated  with  the  much  weaker  wave 
disturbances  generated  during  the  earlier  period,  t  =  0(1),  so  that  the  initial  value  of  all  variables 
( p,p,T,u )  is  zero.  In  contrast,  the  fluid  expelled  from  the  hot  spot  boundary  observed  in  Figs.  5 
and  6  provide  the  boundary  condition  that  drives  the  disturbances.  In  conclusion,  there  will  be  a 
sequence  of  acoustic  disturbances  of  increasing  amplitude  propagating  through  the  unbound  gas 
beyond  the  hot  spot.  Further  solution  development  for  the  acoustic  problem  is  anticipated. 


3.2  Compressible  Limit 

The  rescaled  variables  appropriate  in  the  compressible  limit  are  also  implied  by  the  temporal 
nonuniformities  in  (34)-(36)  when  sTA2  =  0(1). 


(  1  3 

r  _  Tat  ,  Ta-0 

\£  J 

(49a, b) 

(p.T)  =  TA(p.T) 

(49c) 

p=p 

(49d) 

u  =  sTau  =  O  (ta2  ) 

(49e) 

It  follows  that  the  conservation  of  mass, 

momentum,  energy  and  species,  as  well  as  the  state 

equation  in  (l)-(4)  can  be  written  as 

Pi  +  V  •  pu 

(50) 
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(51) 


p(uf  +  u  •  Vu)  =  — —  +  j3T"+n)V 

^(fi+i.vr)  =  -,w.i+rAH>e+/?rr>[^|:+^ 

Yf+u-VY  =  -TASYe~yf  +  pT™  (p/ pSc\ 


p  =  pT 


Q  =  Qe(r,0,(p,t ) 


A HplYe-vT 


(52) 

(53) 

(54) 

(55) 


AH 

where  the  transport  properties  are  identical  to  those  below  (42)  and  8  =  0{1/Ta),  - —8  =  0(1). 

A  He 

The  internal  Mach  number, 
sT,2 1  u 


m=^=^=5te=°K2)=o(  d 


T\'24f 


(56) 


Equations  (50)-(55)  describe  a  fully  compressible,  reactive  flow  with  primarily  small 
transport  effects  except  where  velocity,  temperature,  and  concentration  gradients  are 
exceptionally  large.  The  internal  Mach  number  is  substantial  and  can  be  supersonic.  The 
internal  dynamics  of  the  reactive  hot  spot  must  be  determined  by  a  numerical  solution  of  (50)- 
(55).  The  consequences  of  those  dynamics  will  determine  the  gasdynamic  disturbances  in  the 
far-field.  Kassoy  [2]  has  shown  that  the  inert  analogue  of  (50)-(55)  is  compatible  with  a  strong 
blast  wave  in  the  unheated  external  environment.  In  fact,  it  can  be  demonstrated,  using  the 
results  in  [2]  that  the  blast  wave  is  actually  initiated  within  the  inert  near-field  volume,  far  from 
the  origin  and  at  a  large  value  of  the  relevant  time  scale.  One  can  expect  similar  properties  for 
the  near-field  gas  dynamics  in  a  reactive  hot  spot.  Numerical  solutions  are  required  to 
understand  the  details  and  will  be  the  subject  of  a  future  communication. 


33 


4.  Context,  Summary  and  Conclusions 

4.1  Context 

Transient,  spatially  distributed  combustion  in  a  turbulent  flow  is  thought  to  be  the  source  of 
acoustic  and  stronger  mechanical  disturbances  in  a  LRE  chamber.  Surprisingly,  a  quantitative 
model  for  mechanical  wave  generation  in  a  transient,  spatially  distributed  reacting  flow  does  not 
appear  to  be  available  in  the  technical  literature.  Kassoy’s  thermomechanics-based  formulation 
[2]  for  related  processes  in  an  inert  gas  identifies  a  quantitative  relationship  between  time 
resolved,  localized  power  deposition  into  a  finite  volume  of  gas  and  the  magnitude  of  acoustic  or 
stronger  disturbances  generated  in  the  neighboring  unheated,  unconfined  gas.  The  results 
demonstrate  that  the  thermomechanical  response  of  the  heated  gas  depends  on  a  complex 
relationship  between  the  magnitude  of  energy  deposited  and  the  time  scale  for  that  deposition, 
but  not  simply  on  the  power  deposition  alone.  The  theory  is  limited  to  heating  time  scales  small 
compared  to  the  acoustic  time  of  the  volume,  a  condition  necessary  for  the  pressure  to  rise  with 
temperature  (while  the  density  change  is  marginal).  When  the  energy  deposition  is  quantitatively 
limited,  nearly  constant  volume  heating  occurs  (near-inertial  confinement),  characterized  by  a 
small  internal  gas  expansion  Mach  number  defined  with  respect  to  the  high  temperature  hot  spot 
speed  of  sound.  Temperature  and  pressure  rise  together  within  the  hot  spot.  Localized  high 
pressure,  relative  to  that  in  the  cold  neighboring  gas,  is  necessary  to  create  the  pressure  gradient 
source  of  induced  mechanical  motion  (both  fluid  speed  and  acoustic  disturbances).  Gas  expelled 
from  the  boundary  of  the  hot  spot  (the  “piston  effect”)  is  the  source  of  mechanical  disturbances 
in  the  unheated  environmental  gas.  The  expelled  gas  Mach  number,  defined  with  respect  to  the 
cold  gas  speed  of  sound  will  be  exceptionally  small  when  the  energy  deposition  is  sufficiently 
small  relative  to  the  quantitative  limitation  referred  to  above.  The  resulting  mechanical 
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disturbances  in  the  cold  environment  are  linear  acoustic  waves.  Larger  energy  deposition  within 
the  quantitative  limitation  is  associated  with  much  larger  expelled  gas  Mach  numbers  and 
significantly  stronger  shock  wave  propagation  in  the  cold  environmental  gas.  Beyond  the 
quantitative  limitation  referred  to  above,  the  heating  process  is  fully  compressible,  characterized 
by  a  very  large  internal  Mach  number.  The  asymptotic  theory  demonstrates  that  this  type  of 
thermomechanical  response  to  energy  deposition  is  the  source  of  strong  blast  waves  [23]  in  the 
unheated  external  environment. 

4.2  Summary  and  Conclusions 

The  theoretical  framework  in  Ref.  2  has  been  used  in  the  present  work  to  quantify  the 
thermomechanical  response  of  a  reactive  gas  affected  by  a  localized  exothermic  chemical 
reaction.  The  extension  is  formulated  for  a  subcritical  perfect  gas  with  one-step,  high  activation 
energy  exothermic  kinetics.  The  asymptotic  formulation  demonstrates  that  the  thermomechanical 
response  of  a  reactive  gas  to  localized,  rapid4  chemical  heat  addition  is  similar  in  nature  to  that 
for  an  inert  gas  with  external  heating.  Nearly  inertially-confined  heating  occurs  only  when  the 
chemical  heat  release  and  the  high  activation  energy  are  quantitatively  restricted.  Within  that 
limitation  one  finds  linear  acoustic  wave  generation  for  the  smallest  range  of  energy  addition  and 
stronger  wave  generation  for  more  significant  energy  deposition. 

Numerical  solutions  to  the  equations  describing  the  nearly-inertially  confined  reaction¬ 
generated  heat  addition  process  show  that  a  spatially  distributed  reaction  wave  appears 
spontaneously  in  the  hottest  portion  of  the  hot  spot,  and  propagates  through  the  relatively  slowly 
moving  fluid  at  a  supersonic  speed,  relative  to  the  hot  gas  speed  of  sound.  As  the  front  nears  the 
much  colder  boundary  of  the  hot  spot  it  decelerates  significantly  and  steepens  considerably.  The 
configuration  evolves  toward  a  discontinuous  front  separating  hot,  high-pressure,  burned  gas  on 
4  Heating  time-scale  short  compared  to  the  local  acoustic  time. 
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one  side  from  cold,  low  pressure  reactant  on  the  other  side.  Although  the  pressure,  temperature 
and  reactant  concentration  jumps  across  the  front  are  reminiscent  of  a  reactive  shock  wave,  the 
front  does  not  propagate  at  a  supersonic  speed  relative  to  the  unbumed  cold  reactant  near  the 
boundary  of  the  hot  spot.  In  fact,  during  the  relevant  heating  time  scale  the  entire  process  occurs 
in  a  nearly  incompressible  medium.  The  heating  process  on  the  time  scale  of  interest  has  created 
an  ephemeral,  isolated,  burned  out,  hot,  high  pressure  spot  embedded  in  a  cold  unreacted,  lower 
pressure  gas.  The  large  pressure  gradient  at  the  front  induces  a  local  positive  radial  fluid 
velocity.  Fluid  expelled  from  the  hot  spot  boundary  acts  as  source  of  mechanical  disturbances 
propagating  into  the  neighboring  cold  gas.  The  amplitude  of  those  disturbances  depends  upon 
the  energy  addition  level  during  the  reactive  phase  of  the  hot  spot.  Of  course  there  must  be  a 
subsequent  longer  time-scale  process  in  which  the  pressure  differential  is  relaxed.  One  can 
anticipate  that  a  combustion  wave  will  be  driven  into  the  neighboring  reactant. 

This  type  of  asymptotic  formulation  and  analysis  can  provide  fundamental  insights  into  the 
sources  of  mechanical  disturbances  within  a  LRE  chamber  containing  a  turbulent  reacting  flow, 
consisting  of  spatially  distributed  regions  of  transient  combustion.  In  addition  the  results  define 
crucial  length  and  time  scales,  as  well  as  dominant  physico-chemical  mechanisms.  These 
outcomes  can  inform  the  numerical  modeler  seeking  to  develop  code  capable  of  resolving 
transient,  compressible,  reactive  turbulent  flows. 
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